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FOREWORD 


This study was carried out for NASA by Grumman’ s Advanced Development Office 
under the Shuttle Research and Technology Project (W. Ludwig, Manager; Dr. G. 
DaForno, Aerothermo Mgr.). The contract number is NAS1-11818, March 1972 to 
October 1973. Many of the results presented in this report were initially obtained un- 
der an in-house study. 

M. Rossi developed the numerical method and the DETRAD subroutine package. 
Besides the authors’ efforts, contributions were also provided by A. Jameson, who ini- 
tially suggested the numerical approach and was consultant to the study, and by G. Da- 
Forno, who developed some points on the error amplification and the ‘yes-no’ chart, 
obtained the computer time minimization data, and worked out the extension to temper- 
ature-dependent properties and radiation. G. DaForno also assembled the results and 
wrote portions of the report. 


LATERAL CONDUCTION EFFECTS ON 
HEAT-TRANSFER DATA OBTAINED WITH 
THE PHASE-CHANGE PAINT TECHNIQUE 
By George Maise, Michael J. Rossi 
Grumman Aerospace Corporation 

SUMMARY 

A computerized tool “CAPE” (Conduction Analysis Program using Eigenvalues) has 
been developed to account for lateral heat conduction in wind tunnel models in the data 
reduction of the phase-change paint technique. The tool also accounts for the effects of 
finite thickness (thin wings) and surface curvature. A special reduction procedure using 
just one time of melt is also possible on leading edges. A novel iterative numerical 
scheme was used with discretized spatial coordinates but analytic integration in time to 
solve the inverse conduction problem involved in the data reduction. 

A “yes-no” chart is provided which tells the test engineer when various corrections 
are large enough so that CAPE should be used. 

The accuracy of the phase-change paint technique in the presence of finite thickness 
and lateral conduction is also investigated. 
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NOMENCLATURE 


Symbols 


a 

a U 

A 


c p 

u 


h 

h 

H 

k 

k 

1 

M 

N 

P 

Pr 

q 

Q 

r 

R 

s 

t 

T 

V 

V 

AV 

Ax 

a 

c 

X 

A 

P 

P 


half width of “top-hat” h distribution 
elements of coefficient matrix A 
coefficient matrix 

conduction area between elements i and j 
specific heat 

departure of computed surface temperature at t„, from phase-change 
temperature, Eqn. (11) 

influence coefficient matrix, Eqn. (12) and (14) 

heat transfer coefficient 

modified heat transfer coefficient, Eqn. (2) 

enthalpy 

thermal conductivity; also number of surface elements 

modified thermal conductivity Eqn. (3) 

slab half thickness 

Mach number 

number of elements 

pressure 

Prandtl number 

heat flux 

correction due to TD properties and e ^ 0 

recovery factor 

radius 

distance along surface 
time 

temperature 

velocity along body surface 
injection velocity 
matrix of eigenvectors of A 
volume of element 

distance between elements (center to center) 
angle of attack; also thermal diffusivity 
emissivity 
eigenvalue 

matrix of eigenvalues, also sweep angle 

viscosity 

density 

otjjj/l 2 , nondimensional melting time 
wedge angle 
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Symbols (cont. ) 
Subscripts 


AW 

bg 

CP 

e 

E 

1 

init 

m 

r 

s 

sp 

TD 

(o) 

(n) 

0 

1 

2 


= adiabatic wall 
= background 

= constant material properties 
= edge 
= effective 
= i-th element 
= initial, t = 0 
= melting 

= row number of matrix 
= column number at matrix 
= stagnation point 
= temperature dependent 

= starting value (one-dimensional, semi-infinite) 
= n-th trial of h t 
= stagnation point 
= surface 1 of slab-like geometry 
= surface 2 of slab-like geometry 
= free stream 


Notation 


tV u /V 21 ] = 


Vu 

v 21 _ 


x 



INTRODUCTION 


Subject 

The subject of this study is the data reduction of the phase change paint technique, 
a well known technique (ref. 1 & 2) for obtaining the heat transfer coefficient from wind 
tunnel model tests. Briefly, the model of interest is covered with a paint that melts at 
a known temperature and then at each point on the model surface the melting time is 
measured, i. e. , the time elapsed from the instant the heat transfer coefficient becomes 
steady. The heat transfer coefficients h are then deduced via a calculation of the model 
temperatures with a step-like time history for the h’s. 

Specifically, the situation studied here is that in which the spatial variations in the 
heat loads imposed on the model are so large that lateral conduction in the model must 
be accounted for to obtain accurate heat transfer coefficients. Another situation of the 
same type occurs when lateral conduction is caused by the model geometry; for example, 
lateral conduction is generated in a curved slab of variable thickness under spacially 
constant heat loads. In both situations we place a restriction that it is still possible to 
isolate a portion of the model where the problem is two dimensional. Two such cases 
that occur commonly are (1) a slab (‘thin’ or ‘thick’) impinged upon by a shock, where 
in the direction locally normal to the shock trace the problem can be considered two 
dimensional; and (2) a leading edge (1. e. ) of a model wing or fin, where along cuts nor- 
mal to the l.e. conduction is also essentially two-dimensional. Besides these two com- 
mon occurrences, an arbitrary two-dimensional problem is also considered in this study 
as indicated in fig. 1. 

Two types of data reduction problems can be posed for the situation in fig. 1: 

1) the distribution of the times of melt is given on the surface subjected to heat 
loads and the corresponding h’s are to be derived — this is the typical problem 
in ref. 1 to 4. (Naturally the surfaces over which no data are given are to be 
taken as adiabatic, or if at infinity, at the initial temperature throughout the 
transient.) 

2) the time of melt is given at one point on a leading edge and the entire h distri- 
bution is to be derived — this is a somewhat novel problem proposed by R. A. 
Jones of the NASA Langley Research Center and suggested by the fact that on 
leading edges of very small radius (0. 05 inches for the wing of a typical 1 ft. 
shuttle orbiter model), it is too difficult to determine experimentally the prop- 
agation of the melting line. At most, the minimum time of melt around the 

1. e. can be determined. Approximately, this minimum melt time should occur 
near the maximum heating point or near the stagnation point. Naturally, in 
order to reduce phase-change paint data giving only t,,, at the stagnation point, 
the following information is needed: 

a) stagnation point location in the l.e. section considered, if a* 0. 

b) pressure distribution around the l.e. 
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FIGURE 1. GEOMETRIES CONSIDERED IN THE STUDY 


c) body entropy 

d) T aw distribution around the I. e. , and naturally 

e) a formula for the heat transfer coefficient distribution. 

This information has to be obtained theoretically. To be realistic it is also 
imperative to use: 

a) relatively simple formulae, preferably not involving table look-up or iter- 
ations. 

b) infinite-cylinder approximation for the 1. e. boundary layer. 

Indeed, approximate formulae do exist that appear adequate for the task. How- 
ever, this data reduction procedure has not yet been evaluated, as before this 
study there was no tool for carrying out the calculations. In what follows, this 
data reduction problem is referred to as ‘special problem for the l.e. ’ 

In both types of problems, the normal one and the special problem for the l.e. , it 
is natural to include radiation from the model surface and model material properties 
variable with temperature. These two effects are grouped together because they both 
make the reduction non linear and therefore are handled in a straightforward manner 
only by numerical methods. For the current model materials (Stycast, other epoxies, 
etc.), both these effects are small. In particular the effects of temperature-dependent 
(TD) have been found to be very small both in one-dimensional and two-dimensional 
cases. Note that the TD solution should be compared to a constant property (CP) solu- 
tion where the CP properties are properly chosen. The rational way is to evaluate k 
and c p as the average of the values at T inlt and T m , as indicated in fig. 2. Then, during 
the time history from T lnlt and T m the properties are correct on the average. Of course, 
at later times, when T>T m , the properties are locally less accurate, but to have an ef- 
fect, this inaccuracy would have to influence regions far away where the paint melts 
later. It is this influence that is minimal. For example, in a typical one-dimensional 
case with a Stycast model there was only 0.4% difference between h TO and h cp . Similar 
results were found in two dimensional cases (calculations were done in direct problems 
using AG TAP (ref. 5). In spite of the fact that they are small, that the k and c p varia- 
tions with temperature are not always available and the fact that the properties in- 
homogeneity can be more important, radiation and TD properties are here considered 
for completeness. 


Problem 

Two specific problems are attacked in this study, i.e. : 

i) the development of a numerical method and relative computer code for reducing 
phase change paint data for the three geometries of fig. 1 and for both types of 
experimental inputs mentioned above; 
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ii) the preparation of a ‘yes-no’ correction chart, that quickly tells the test engi- 
neer when the lateral conduction corrections are large enough (say 10% on h) to 
make reduction by ‘semi-infinite slab’ (ref. 1) or finite slab (ref. 3, 4) too in- 
accurate and therefore the use of the computer code above necessary. 

The chart should, of course, be simple. Moreover, the computerized data reduc- 
tion should allow for variable material properties and radiation. Naturally, the code 
itself should be easy and quick to use, especially as far as geometry input (e.g. , grid 
lay-out, elements selection, conduction paths calculations). 

Target Run Times 

To produce a practical computer tool for the data reduction in the presence of (two- 
dimensional) lateral conduction, we set a maximum run time on the CDC 6600 on the 
order of 6 minutes for a typical slab case with h obtained to an accuracy around some 
1 to 2% (which is an adequate value to impose on the numerical method, as it is an order 
of magnitude smaller than the absolute accuracy of the experimental h). 

Prospects for Data Reduction in 
Three-Dimensional Geometries 

Currently, the maximum service temperature of model materials such as Stycast 
imposes various limitations on the test conditions the model can be tested at, the max- 
imum run time of the test and the number of times a model can be used. Metal, e. g. , 
stainless steel models would eliminate these drawbacks of the epoxy models, at the 
price (among other things) of a much larger conductivity and therefore fully three-di- 
mensional lateral conduction. Reducing data in the presence of 3D lateral conduction 
effects would require in the numerical method a large number of elements to describe 
the model (or a 3D portion of it). While not a part of this study, such a prospective 
extension of the numerical method should be kept in mind. As a target, one may per- 
haps put the maximum run time, on the CDC 6600, at 30’ . Of course, the feasibility 
and practicality of such an elaborated data reduction would have to be judged on many 
more counts than the mere low computer time. 
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NUMERICAL METHOD FOR THE LATERAL CONDUCTION 


INVERSE PROBLEM 


Possible Approaches 

The most obvious approaches that could be taken, are finite difference methods 
(implicit or explicit) or finite element methods for the space variables and time, or a 
mixture of finite element for the space variable and finite differences for time. These 
are the methods used in the current computer tools for direct heat conduction problems, 
typically SINDA (ref. 6), AGTAP (ref. 5), etc. Along these lines the one-dimensional 
inverse problem has been also solved successfully (ref. 4.) 

Unfortunately, none of these conventional approaches (finite difference or finite ele- 
ment, implicit or explicit) seems to offer much hope for a practical tool for two-dimen- 
sional inverse problems because of the extremely large machine times required. The 
basic reason is that, in iterating on the h^ once a set of h t is guessed, one needs influ- 
ence coefficients of the type 8T t /8hj and to obtain these coefficients it is necessary to 
calculate a very large number of direct problems. To try to put the problem on quanti- 
tative, even if approximate, terms, figure 3 shows some comparisons of running times 
for explicit and implicit methods coupled with numerically determined influence coeffi- 
cients (Newton’s method). For the purpose of comparison, it was assumed that all 
methods converge in three iterations , which tends to put the method developed below at 
a disadvantage. * 

The explicit method suffers from the fact that for physically small elements the 
critical time step beomces very small. Furthermore, if stainless steel models are to 
be considered, the critical time step is further reduced by a factor of 12 over that of 
Stycast with the corresponding increase in computer time. Thus, the interesting con- 
clusion follows, from fig. 3, that the explicit method has to be ruled out. 


♦The curves in fig. 3 were determined as follows. For the explicit method, we ran a 
121 -element Stycast problem (heated on one side) in the direct mode with a thermal 
analyzer (ref. 5) . To do the corresponding inverse problem with the thermal analyzer, 
the influence coefficients would be determined by perturbing all the h t ’s in turn and re- 
running the program. This would lead to a total running times of 21 minutes. The in- 
dicated proportionality of computer time with n 2 (n number of elements) would result if 
the number of surface elements to the total number of elements were kept in the ratio 
of 1 to 11. For the implicit method we timed one of our very efficient finite-element 
structural mechanics programs which resembles (mathematically) the implicit heat 
transfer problem. To do the corresponding inverse problem the same arguments that 
were used for the explicit method apply. A variation of computer time with n 2 was as- 
sumed because the implicit method involves matrix operations. Finally, for the meth- 
od developed here, fig. 3 presents running times actually clocked for an early, non- 
optimized version of both method and code. 
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Computer time with the implicit method depends directly on the time step that is 
chosen. Figure 3 assumes that for accurate computations 50 time steps are required. 
The computer time could, however, be lowered at the sacrifice of accuracy in the time 
integration. However, the basic method does not possess the potential for significant 
reductions in computer time, and therefore, it is very doubtful that it could ever handle 
three-dimensional geometries within the 30' target time. 

Method Developed Here 

In order to deal with the difficulties just mentioned which are characteristic of the 
problem at hand, a different idea (suggested in this context by A. Jameson) was devel- 
oped as follows. 

Imagine to discretize as usual the structure, in lumps for example in the spirit of 
the finite element methods. Write the ordinary differential equation governing the tem- 
perature history in each element. The idea is now not to discretize the time variable 
and carry out the integration in time numerically, but instead to leave this variable con- 
tinuous and carry out the integration analytically. The method therefore belongs to the 
class of the semi-descretized variables or the hybrid analytical-numerical methods. 
Hybrid methods have not yet received much attention, even though in the literature a 
few cases have been reported. Their potential seems to lie in their ability to^ incorpo- 
rate the best aspects of both numerical and analytical methods. In this essential aspect 
the method developed here differs from the usual finite-element methods used in heat 
conduction. 

The great advantage in not discretizing the time emerges when the inverse problems 
of phase change paint data are to be solved. These inverse problems are described by 
multipoint boundary value problems in time. Since the integration here is done analyti- 
cally, the result is an explicit expression for the temperatures of all elements as func- 
tions of time and all the heat transfer coefficients, i.e. : 

Tj = T t (t, hj, h 2 . . . h k ) 

When the t^ and T m are given, this equation has to be solved iteratively for the h^s. 

An efficient iteration requires the influence coefficients G tj = 8T t / 9hj. Two approaches 
are avilable to obtain these influence coefficients. First, since there is an analytical 
expression for T t , the G tj can be obtained directly by differentiating this expression 
with respect to all the h t ’s. With any numerical scheme that discretizes time, the time 
integration must be repeated k times (k = number of surface elements) to obtain these in- 
fluence coefficients. In other words, the usual methods, whether implicit or explicit, 
require k+1 integrations in time, while our scheme needs just one. As the number of 
points increases and also the melt times become larger (driven by the quest for higher 
accuracy in the data reduction, larger model conductivity and higher model service tem- 
peratures), k integrations in time are bound to cost, in terms of computer time andstor- 
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age, more than the calculation of an explicit expression for G^. The other approach is 
to calculate the G tj , from the exact ‘analytical’ or ‘infinitesimal’ expression, from the 
approximate finite-difference expression: 

_dTi ,_ T 1 (h 1 ,h 2> h 1 + Ah 1 ....h k )-Ti(h 1 h 2 ...h k ) 
ij 8h j Ahj 

with a small but finite Ahj. This requires one extra set of eigenvalue calculations per 
surface point. Therefore in this case we trade the k+ 1 time history calculations of the 
conventional methods with k + 1 eigenvalue calculations (k of which have only one h per- 
turbed) . 

In both cases, then, this unconventional analytical-numerical method permits the 
substitution of the numerical integration in time by one or more eigenvalue calculations 
and opens the possibility of taking advantage of the considerable amount of work done 
recently on fast eigenvalue calculations in the areas of aircraft control theory and struc- 
tural analysis*. This gives latitude for considerable advances in computer time, stor- 
age and size of problems handled. Such latitude is not evident in the typical heat con- 
duction methods . 

Naturally, all this needs to be put on a quantitative basis. Figure 3 anticipates 
some indications of the machine time requirements for the method developed here. One 
point to be kept in mind is that these times are for a non-optimized version of method 
and code and that the method itself has considerable potential for time reduction. 

Clearly, this hybrid analytical -numerical method can be carried out within the 
framework of finite element approaches or within the framework of finite difference ap- 
proaches as far as the space variables are concerned. For the problem at hand, there, 
may be only minor differences between these two approaches in the final discretized 
equations. We patterned our space discretization after the simplest version of the finite 
element approach, namely uniform distribution of each quantity within each element, in 
the belief that what is most important for large improvements in the numerics is an 
imaginative approach on the broad issues rather than sophistication in detailed matters. 
However, one point of the spatial discretization has been paid attention to, namely 
boundary elements are treated in a (very simple) way that assures — within the uniform 
distribution assumption — maximum accuracy in imposing the boundary conditions, as 
will be seen below. 


*The method developed here was transplanted from control theory, but most of the rou- 
tines and eigenvalue-eigenvectors numerics originate from structural analysis. 
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Equations for the Basic Case 


The key equations in our method are now presented. We give in this section the 
basic case of variable T AW , no radiation, times of melt that are all given. 

After the model’s structure, or more likely, a portion of it is suitably subdivided 
into n elements, the temperature response of each element i can be represented by 

pc p Av * 7k = hiAi {Tm > 1 " + S k ^ (T J “ Ti) (1) 

(hi term suppressed for interior elements; j elements adjacent to i share with i the in- 
finite area A^). The symbols are defined in the nomenclature. For simplicity, re- 
define new coefficients. 

hi = h^ (2) 



The equation can be rearranged to read 


pC p AVl ~dT = ~ (hi + & k ‘J )T ‘ + § k ii T l +h i T AW,i (4) 

These equations constitute a system of n first-order linear differential equations with 
constant coefficients. Because of the variable T AW , these equations are inhomogeneous. 
Writing (4) in vector notation 

=BT + F, T(0) = T lnlt 
or alternatively 


^ = M- 1/2 AM 1/2 T + M _1 F f (5) 

dt > 

T(0) =T lnit ) 

where M is the diagonal matrix made up with the pC p AV^i, A= M -1/2 BM~ 1/2 is a symmet- 
ric matrix where 


b 

b 


u = - (hj + ^kij), hi = 0 if i is interior element 


( 6 ) 

(7) 
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and F is a constant vector 


F = 


hiTAw, i 
hiT AWi j 

h k T AW>k 


( 8 ) 


The solution to the initial value problem (5) is written in the usual fashion as 

T = T 0O + M' I/Z Ve At V T M 1/2 (T lnit - T J (9) 

Where the t — 00 temperature T„ is the solution of the RHS of (5) set equal to zero: 

M- 1/2 AM 1/2 T oo +M j F = 0 (9a) 

and where A is a diagonal matrix formed with the eigenvalues of the matrix A, and V 
is the matrix of eigenvectors of matrix A, V T is transposed of V. f Thus, once the ei- 
genvalues and eigenvectors are determined with one of the standard subroutines, the 
temperature of any element at any time t is readily evaluated by eq. (9). For the prob- 
lem at hand, only the temperatures of surface elements at t^j are of interest. 

Equation (9) solves the direct problem, since it gives the temperature distribution 
when T init and hj are given. The problem of interest is really one of determining the 
parameters h t from a partial knowledge of the direct solution (T = T m>I at t = C m 1 1 for each 
surface element). Therefore, the solution of the inverse problem is given by equations 
(9a) and (10). 

T m,i = 6i'M- 1/2 Ve Atm ' 1 V T M 1/2 (T lnit -TJ+T 0O ) 

A, V and T« are functions of hj ) (10) 

Here 6? is the i-th row of the unit matrix of n elements. Equation (10) is a system of 
transcendental equations in the unknown iq. By solving it one determines hj. 

In order to solve eq. (10), an iteration in the hj is used. As starting values for the 
iteration, the hj are first computed using the one-dimensional semi-infinite slab theory 
(ref. 7). The temperature of the surface elements corresponding to the given t^j are 
computed by eq. (9). This is of course the RHS of eq. (10). These temperatures should 
all equal T m . Considering h a vector parameter, the errors in the temperature are thus 

fi® = T m -6fT(t m , i ) (11) 

To obtain new values for the hj’s, use is made of Newton’s method of iteration. It 
is necessary to evaluate the influence coefficients G^, (G in the matrix form), where 

i f n 

t If T AW is constant in space, T* = T AW is also constant. If H = 0, T*,^ (2 m k T inltfk )/ 

If n k 

(Zm k ) is constant which of course, for T lnit constant in space, becomes T lnlt the same 
constant. These results show that the discrete model is consistent. 
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( 12 ) 


iy. 


9f. 

Gl J _ "ah " 6l 


. sT iy 


ah, 


Since the T(t m>i ) are known, the derivatives, 
The result is 

G u = 6 ? M 1/2 Vtmj 1 Qij(t mfi ) V T M 1 / 2 (T lnit - TJ + 

e^m , i — e*s 4 m , i 


Qjrs~ 


\V jr V JS ' 






X r ^ Xg 
Xr — X3 


i 


61 = [ 0 *— 0 l 0 »«* 0 ] 


A„ = 


i 

o *»* *0 ••••0 
o* # oio***o 
()•••• 0 ••••0 


i.e. , the G u can be calculated formal- 
( e Atm,i _ DA- 1 V T M- 1 / 2 A jj (T AW -TJ] (13) 


Although the expression for G tJ appears to be quite complicated, it is to be noted that 
V and \ i are already available from the calculation of the T(t m i ) (the RHS of eq. (10)), 
and the multiplication of terms in eq. (13) can be performed in such a way that matrix- 
by-matrix multiplication is avoided. 


There is another possibility for calculating the G 1} . Rather than the infinitesimal 
expression (13), one could use the approximate finite-difference expression: 


_ 9Tj Tj (hj, h2. . . . hj + Ahj. . . , hjj) Tj^ij, I 12 h^) . . 

o u -_- — U41 

ohj Ahj 

with a small but finite Ahj. This requires the evaluation of one extra set of eigenvalues 
and eigenvectors (for the perturbed sets of h } ) for each i and the corresponding Tj at the 
time of melt t m>i . Naturally there is no point in using a central difference that will re- 
quire two extra sets of eigenvalues and eigenvectors. Later on, we will discuss the 
question of which of the expression (13) or (14) for the G u is most convenient. 

Once, however, the G^ are determined during each hj in the following way. 

For a change 6 h in the parameter vector, eq. (12) gives 

6 f = G 6 h + 0 II 6 h II 2 (15) 

To make the errors f vanish at the next trial f should equal — f from the previous trial. 
Ignoring the higher order terms, eq. (15) can be solved for the required change in K: 

6 H=-G- 1 f (16) 
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To carry out this step requires an inversion of the G matrix. It is to be noted that this 
matrix is typically much smaller (it is equal to the number of surface elements on which 
h is unknown) than the A matrix. The updated set of hj can be formed as 

6h J <1) = h/ 0) + 6h j (17) 

The above procedure is repeated until eq. (10) is satisfied. When this occurs, the de- 
sired hj’s have been determined. 

As will be shown below, the iteration scheme just described is very powerful. For 
example, convergence to 0. 1% error is generally achieved after three to four iterations, 
for the most severe case, when the times of melt at all surface points are given. 

Temperature Dependent Properties and Radiation 

Radiation and TD properties* can be simply included in the numerical method that 
we have developed. This is not quite obvious at first sight since the method depends 
upon the analytical solution for the ordinary differential equations and we are able to 
write the analytical solution only for problems linear in time. The usual radiation for- 
mulae or material properties variable with temperature make the problem non-linear. 

The analysis can be extended to handle non-linearity in a manner that preserves the 
power of the numerical method. The extension amounts to “lag behind" in the non-lin- 
ear terms so as to build up the exact solution through a succession of iterates, each of 
which is linear. The iterates are not new ones, rather the ones already necessary for 
h 3 . The key point is that the analytical time integration for each element is maintained. 

With radiation and variable material properties, the starting equation replacing 
(1) is: 

pC vt (T t ) AV t hjAj <T AW|t - T t ) - aei - T 4 B G , i )A l 

+ .Ek ij Cr i ,T j )^(r j -T i ) (is) 

with the same convention as in (1) to suppress the first two terms in the RHS for interior 
elements and the same conventions for j. 


^Naturally, radiation with arbitrary emissivity and background temperature for each 
surface; and also arbitrary C p (T) and k (T), supplied, for example, in the usual tabu- 
lar form, p is taken as constant. 
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Using constant but suitable values for c p and k 

- = c p (T lnit) + c p< T m) 
p 2 


jT- k(T tBlt ) + k(T m ) 

2 

the equation is now recast as 


HT — A, , 

p ^ v ‘d r =h » A ‘ (T ^.i - T t> + ^ k (t > - T * ) 

uu j5t j 

+ Q t (T,(t)jTj (t)) (19) 



( 20 ) 


Naturally Q t summarizes all the non-linearity in time and is made up of three effects as 
is clear from inspection of (20). 


Proceeding as before, equation (5), (9), (10) and (13) are now replaced by (consid- 
ering Q a known function and maintaining the same definition of T*): 

dT 

— =M* 1/2 AM 1/2 T + M- 1 F+Q (6') 

T=M- 1/2 Ve At V T M 1/2 (T lnlt - Tj + f* M‘ 1/2 Ve A(UT) V t M 1/2 Q(t) dT + T. (9') 


T M ,i =6^M- 1/2 Ve Atm *i V T M 1/2 (T lnlt - Tj 

+ M- 1/2 Ve A (t m,i-T) v T M 1/2 Q t (T)dr + T*, 

G ij =^[m- 1/2 VQ j (t^) V T M 1/2 (T lnlt -Tj 
+M" 1/2 VQj (t m>1 ) V T M 1/2 Q t (^ tl ) 

1 

Ohj J 

(h m = 0 r*s 

matrix H < - . 

1 e x r tm »i_l 


(10') 


(13') 
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In these equations all the first terms are clearly those obtained in the absence of radia- 
tion and material properties variation. In particular, the equation solved for the in- 
verse problem looks simply 

T m ,i = (T mtl ) NO RAD, + / 0 tm,1 M- 1/2 Ve A(t «M^V T M 1,2 Q 1 (T)dT (10") 

CP 

This equation that gives the h £ cannot be solved alone any more, but must be solved to- 
gether with the integral equation (9') that gives the T £ (t). However, recalling that radia- 
tion and TD properties are not dominant effects, Q is relatively small and therefore an 
iteration in Q should rapidly converge. Therefore, in each h £ iteration, Q(t) can be 
taken approximately as the value computed a posteriori from the temperatures of the 
previous h £ iterate. In the first iteration, Q(t) = 0. Of course, as h £ converge to the 
exact values, so do the Q(t)s. 

In this way linearity is maintained in each h £ iterate and the calculation method re- 
mains as indicated in the previous section. No new iteration for Q is needed. The only 
addition is the calculation of Q for equation (10") and for this temperature history up to 
t^ is needed from equation (9'). 

The integral in (10") is discretized in the most straightforward fashion 

T„,i “ (!*><** + l, £ * Y 1 [M ‘ 1/2 Ve A(t m,i-n + i>v T M 1 / 2 Qj(r 1+1 ) 

+ M- 1/2 Ve Attm » 1 - Tl, V T M 1/2 Q i (T 1 )] (21) 

Ar=-^ 

Tnajc ~ 1 

Tj= (1 — 1)At 

lnuut arbitrary 

Equation (20) gives Q t (t £ ) via the T(tj) from equation (9'), which is discretized just as 

( 21 ). 

No great sophistication is needed in the number of time intervals l mstt nor in the 
discretization of the integral near t = 0 since the contribution of Q is altogether small. 

In the same spirit, when using the analytical expression for the G tJ , we can retain 
(13) instead of (13') even with radiation and TD properties. Experience with this meth- 
od has been that only approximate values of the G tj are needed to drive the hj iterates 
to the solution in a few of iterations. Therefore, there is no need to add the small ra- 
diation correction to the G 1£ whose only function is to drive the iteration. 

Numerical results show that under normal conditions, radiation and TD properties 
do not affect the number of iterations and have very small effect on G tj . However, the 
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extra calculations do cause an increase of computer time by a factor of about 2 for radi- 
ation and 3 for radiation plus TD properties. Naturally a method that integrates numer- 
ically in time would entail a small increase due to these two effects. However, since 
the effects appear to give a maximum change of the order of 1% in hj for current situa- 
tions, the present method allows to take advantage, for the majority of the applications, 
of the fastest version of the method. This advantage does not exist in the method re- 
quiring the numerical integration in time. 
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Geometry Discretization and Calculation Of 
Capacitances and Conductances 


In this section we document the way in which the three types of geometries of in- 
terest (fig. 1) are discretized into elements and how capacitances and conductances 
are calculated. The calculations (grid lay-out and capacitances) are done in a stan- 
dard way since they are automatized in the computer code so that the tedious labor of 
input preparation typical of multi-dimensional conduction codes is mostly eliminated. 

The grid lay-out for the slab, the l.e. and the arbitrary four-sided two-dimensional 
geometry is illustrated schematically in figures 4, 5, and 6 respectively. 

Some points to be noted are: 

a) Control points. The temperature points are taken to be in the middle of the 
elements except for the surface elements, where the temperature points are 
taken on the surface. This procedure results in greater accuracy in surface 
temperature calculations, which is particularly important in phase-change 
paint applications. It eliminates the needs of extra-fine resolution near the 
surface. 

b) Slab thermal capacitances and conductances are calculated as follows (see fig. 
4): 

Ci = Pc p Ax 1 Ay 1 (22) 

and 

<23) 

c) Leading Edge. Fig. 5 represents a geometry in a plane normal to the l.e. 

The cylindrical portion of the geometry is divided into elements by concentric 
circles and by rays. The wedge portion is divided into rectangular elements 
except near the centerline where the elements are trapezoidal. Conductances 
such as E-H, etc. , see fig. 5, are taken to be zeros since these elements 
meet at a point. In the nose region conductances between elements are com- 
puted by the usual logarithmic relationships. For example: 

= (24) 

P 1 a 

and 


k A -B = k^ln -1 ^ 

r e 


(25) 
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Ax, Ay ARBITRARY 

NUMBER OF ELEMENTS IN BOTH Dl RECTIONS ARBITRARY 
2 EXAMPLES SHOWN WITH ACTUAL DIMENSIONS 


Ax U 

ARBITRARY j 



0.380” 

(9.65mm) 

(a) "SEMI-INFINITE” SLAB, 121 ELEMENTS 



LATERAL GRADIENTS OF MELT TIMES: 
HIGH Ax i LOW 



FIGURE 4. GRID LAYOUT FOR SLABS 
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Ax, Ar ARBITRARY 
NUMBER OF ELEMENTS (IN THE 
r AND x DIRECTIONS AND ALSO 
IN THE CIRCULAR PORTION 
OF NOSE) ARBITRARY 
0'S EQUAL 



THE OTHER HALF EXACTLY 
THE SAME IRRESPECTIVE 
OF ANGLE OF ATTACK OR 
GRADIENTS IN INPUT 


FIGURE 5. GRID LAYOUT FOR LEADING EDUES 
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One complication arises near the center of the circle. The conductance be- 
tween adjacent elements such as F and G becomes infinite according to the 
above relationship. To circumvent this problem a small “hole” equal to 0. 01 
r CAP i s assumed for the center. The conductances between rectangular ele- 
ments are computed in the same way as for the slab. For the conductances be- 
tween dissimilar elements (e. g. , I-J in fig. 5) are .computed in two parts and 
added according to the conventional rules. 

d) Arbitrary Four-sided Geometry. The geometry contour is assumed to be 
given pointwise. The grid lay-out is also specified pointwise as indicated in 
fig. 6. To compute by hand the conductances between these elements in this 
case would be an extremely tedious and time-consuming task. This task can 
be done automatically, when the coordinates (x,y) of the corners of the ele- 
ments are given. To compute the conductances between arbitrary quadrilater- 
al elements, use is made of the analysis given by Dusinberre (ref. 8). Dusin- 
berre presents the relationships for computing conductances between triangu- 
lar elements. Since all quadrilaterals can be subdivided into triangles, these 
relationships carry over to the problem at hand. Basically, two situations 
exist as shown in fig. 6. If the obtuse angles of the quadrilateral are on oppo- 
site corners, the resolution into acute triangles is accomplished by passing 
a diagonal through these corners (case 1). If the obtuse angles are adjacent, 
the quadrilateral must be divided into three triangles (case 2). The resulting 
conductances are given by the following equations: 

case 1: 


*A + *B 

where 

k A = i k(cotDAC +cotADC +cotABC +cotACB) 
and 

k B =^k(cotCBF +cotBCF +cotBEF +cotBFE) 
case 2: 

-t 

kj * J+1= T7IL 

and 


(26) 


(27) 

(28) 


(30) 


k c = £ kicotBAG +cotAGB +cotCDG +cotDGC +cotBCG +CBG] 


21 



(31) 


k D = 5 k [cotE BH +cotBHE +cotFCH +cotCHF +cotHEF +cotHFE] 

The capacitance of the arbitrary quadrilateral element (for example i+1 in 
the insert of figure 6 ) is given by 

c M = 5pc p [(x G - xc) (y H - y B ) - (x H - x B ) (y G - yc) 1 < 32) 

The code contains the logic to examine the element and decide on the appropriate sub- 
division into triangles. The code then proceeds to compute the capacitances and con- 
ductances as illustrated above. 

Special Problem For The Leading Edge 

As mentioned in the introduction, in this type of reduction problem one obtains 
from the model test just one melt time, the minimum melt time around the 1 . e. , while 
the h distribution is obtained from theory. 

The data reduction problem then is simplified to one of determining only one un- 
known, i. e. , the magnitude of h at the stagnation point h 0 . In performing the heat con- 
duction analysis, the complete geometry must, of course, be treated, however, the 
iteration to convergence is performed on only h 0 . The given time of melt is most con- 
veniently chosen at the stagnation point in the neighborhood of which melting is expected 
to occur first. 

The problem is solved as usual in a plane normal to the leading edge. The infinite- 
cylinder approximation is used. Approximations that are also accepted in order to ob- 
tain simple formulae are: perfect gas with an effective y, constant Prandtl number, 

H<xT and ‘cold walls’. The boundary layer is taken as laminar on the entire 1. e. 

The equation used fork/ho is the Lees’ formula (see for example ref. 9; assump- 
tions are cold wall and m^T): 


_h F 
ho 

where (the meaning of the symbols may be clarified by reference to figure 7): 


P 

P 02 ^ 




(33) 


(34) 


Of course this equation is written in a plane normal to the 1. e. 

To evaluate F quantitatively, one needs the distribution of pressure and velocity 
around the leading edge, and the position of the stagnation point. The pressure dis- 
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tribution is evaluated with the modified Newtonian law (between the stagnation and the 
two sonic points) joined smoothly to a Prandtl-Meyer expansion downstream of the sonic 
points. The modified Newtonian relation for the pressure coefficient is 

Cp = Cpo sin 2 $ (35) 

Expressing C p in terms of pressure, the pressure is given explicitly by 

p = Pooj^-^r — 1^ sin 2 <p + lj (36) 

where p 02 is the pressure behind the normal shock. From (36) the velocity at the edge 
of the boundary layer Ug follows by Bernoulli’s equation: 



These expressions are used as long as Ug remains subsonic (see fig. 7). Beyond the 
sonic points the Prandtl-Meyer relations are used. The wall turning is related to the 
Prandtl-Meyer function in a very simple fashion. For expanding flow one has 


V 2 = Vi + I <p 2 - <Pl I 


(38) 


where v is a unique function of Mach number 


v(M) tan^j^- (M 2 - lT-tan" 1 ^ 2 - : 


(39) 


Thus as one follows around the contour of the leading edge, v is given by eq. (38). 

From v the Mach number is determined by a simple trial and error application on eq. 
(39). If one again assumes isentropic flow, all other flow properties e.g. , p, u are 
readily related to their corresponding values at sonic conditions through the Mach num- 
ber. 


After p and Ug are determined around the body, F is evaluated by numerical inte- 
gration. One can then substitute F into eq. (33), divide by \f (dug/d s) 0 and obtain the 
desired ratio of heat transfer coefficients. It is actually not necessary to evaluate 
(dUg/ds)o. Since h/h 0 must equal 1.0 at the stagnation point, V (dug/d s) 0 is simply equal 
to F at the stagnation point. 

For the stagnation point position, the Newtonian rule is used, namely, that the 
stagnation point is where the “free- stream” velocity (i.e. , the velocity vector in the 
plane normal to the 1. e. ) intersects the 1. e. contour at 90°. 

The distribution of the adiabatic wall temperature is also required. It is deter- 
mined by the simple approximate relationship 
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FIGURE 7. SCHEMATIC OF THE FLOW FIELD CALCULATION AROUNU 
I EADING EDGE FOR THE SPECIAL PROBLEM 




(40) 


VPr = 


Taw~ T e 

T 0 ~T e 


or 


T AW =/P?(T 0 -T e )+T e (41) 

where Pr is given. 

In conclusion, if one wishes to make use of this special reduction procedure for 
l.e. ’s, the additional input information needed consists of the free stream M 0Ofetf# 7 ',p 0O 
T 0 , Pr, the gas constant R and the angle of attack. These quantities are to be under- 
stood given in a plane normal to the l.e. 

The many approximations embodied in the formulae chosen raise the question of 
whether or not these formulae are adequate, even though it is almost unavoidable to 
keep the complexity of the problem within reasonable bounds. It is worthwhile there- 
fore to state briefly how accurate the approximations are: 

• Stagnation point and point of maximum heating can differ considerably so that 
the position of earliest melt need not be the correct position to start the in- 
tegrals in equation 34. Angular displacement can be of the order of 20 to 30° 
on a swept leading edge at higher a (for example, ref. 10). Unfortunately, the 
the simple Newtonian rule for the stagnation point position is not accurate in 
such cases because the effective angle of attack 

tana E = tana/ cos A (42) 

is very high (for example for the wing l.e. of the NASA MSC 040A orbiter con- 
figuration during reentry, a E ~50°). Angular errors of the Newtonian formula 
can easily be 30° to 50° (ref. 10). We are not aware of a more accurate for- 
mula for determining the stagnation point position on leading edges. 

• None of a dozen or so approximate methods for obtaining the pressure distribu- 
tion around the 1. e. is accurate in the shoulder region, including the typical, 
not-too-complex modified Newtonian plus Prandtl- Meyer. Typically C p can be 
in error by 50% in the shoulder region. Pressure gradients cannot be used at 
all. 

• When comparing with exact inviscid calculations (obtained with computer code 
of ref. 11), the entropy on the body is found to be well predicted by normal- 
shock entropy with effective free-stream Mach number. For a planar wing 
with dihedral <p, sweep A, angle of attack a 

M eff = M^cosAe 

sinA E = sinAcosa +sina cosAsin$ 
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Refinements on these items should not be included in the method at this point. The 
approximations involved have to be kept in mind, but a realistic first step is to use for- 
mulae such as (33) and (41) carefully avoiding explicit pressure gradients. This will 
already represent a considerable step permitting the evaluation of this special data re- 
duction method. 


As far as the numerical method used to solve the equations, there are only minor 
modifications to the Gjj equations (eq. (13) and (14)) and to equation (16) for the updates of 
the h^. Since there is only one h 0 , we have instead of (14) 


, 3T 0 ^ T 0 (h() + Ah, hj +Ahj h k + Ahj.) — To(hg, tq...h k ) 

Aho 

Ahj = ^4 Ah 0 (14") 

ho 


Similarly (11) and (16) are replaced by: 


fo ~ I'm ~ To(tm> o) 
6 hj =“= I f o 

hg u oo 


Minimization of Computer Time 


( 11 ") 

(16') 


Areas for Computer Time Minimization 

In writing the equations of the method, some alternatives were left open in the way 
the calculation is carried out. The machine time required by each alternative will now 
be investigated and the most attractive alternative selected. Specifically, the areas in 
which there is still some freedom left are: 


a) whether to calculate all the eigenvalues of the A matrix (see for example equa- 
tion (10)) as opposed to calculate only the dominant ones and therefore save 
time in the eigenvalue-eigenvector (E&E) calculation; 

b) whether to calculate all the terms in an equation such as (10) or limit the cal- 
culations to only the eV terms containing a dominant eigenvalue; here time is 
saved in the matrix operation to calculate the temperatures and the G^; 

c) what method to use for calculating the E&E’s; 

d) finally, how to calculate the Gjj, whether analytically through the explicit for- 
mula (13) or numerically through (14). 

A few comments will illustrate these items. 


As far as items a) and b) are concerned, this matrix reduction technique has been 
used successfully in vibration analyses (e. g. , ref. 12, 13 and 14) where problems in- 
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volving several thousand elements are common. When applied to heat conduction, the 
idea amounts to representing the transient thermal behavior of the interior elements 
by a smaller number of “generalized” coordinates. For example, the temperature of 
each element i is normally written as the sum 

b n t. * 

T i (t)=EPi j e X 3 t 

The number of eigenvalues and, consequently, the number of terms to be added is exact- 
ly equal to the number of elements. But many of the terms in the sum are so small that 
they could be neglected without loss of accuracy in T t . Table I gives for a typical case 
the error found by keeping only the largest terms. Clearly, sufficiently accurate values 
of Tj can be obtained by reducing the A matrix even by a factor of 10. The deletion of 
the remaining terms may thus result in a significant savings of computer time, partly 
because the subdominant eigenvalues need not be computed (item a) and partly because 
all the subsequent calculations involving the subdominant eigenvalues need not be carried 
out (item b). One problem, however, to be resolved with item b) is the effect of using 
only the dominant eigenvalues in obtaining the influence coefficients G { ^OTj/Ohj when 
these coefficients are calculated through their analytical expression (13). Because of 
the differential nature of this expression is it not obvious that it will not depend in an 
essential way on the complete set of E&E’s. 

As far as item c) is concerned, standard and efficient methods exist for obtaining 
the E&E’s (all or just a given number of diminant ones) of a real symmetric matrix, 
such as A. But of course, there exists the possibility of using methods that are partic- 
ularly efficient for the special characteristics of the matrix A, which is sparse and 
banded. Another point to be kept in mind in choosing a method is that the eigenvalue 
calculations are repeated for each matrix A, i. e. , each time the matrix A is changed 
during the iterations in the hj. Therefore, iteration methods for the E&E’s of eachmay 
be very attractive because there are already available good starting approximations for 
the E&E’s, those of the previous A iterate. Of course, it is not within the scope of this 
study to embark in the development of new methods for E&E’ s calculations; the aim is 
rather to exploit the best state of the art. 

Finally, as far as item d), the question of how to obtain the G lj( whether through 
(13) or (14), depends not only upon which of the two equations takes up more time, but 
also whether or not equation (13) requires all the E & E’s thereby affecting the best pro- 
cedure selected under items a) and b). 

Items a) and d) were investigated through experimentation via a pilot computer 
code. * The problems taken as typical in the experimentation are: (i) a slab heated 
on two sides with constant adiabatic wall temperature and the special l.e. problem, (ii) 
with constant material properties and negligible radiation, and (iii) with a total number 
of elements of the order of 100. These are the most common current problems to which 
the method is intended and for which minimum computer time should be assured. 


*The code which contained just a few double precision operations was exercised on the 
IBM 370/165. The results were spot-checked on a single precision code on the CDC 
6400 with substantially the same results. 
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TABLE I. TYPICAL ERRORS INCURRED IN NEGLECTING SUBDOMINANT EIGENVALUES 
TERMS IN THE TEMPERATURES CALCULATIONS 


NUMBER OF DOMINANT EIGENVALUES + EIGENVECTORS n e 



2 

3 

4 

5 

10 

20 

30 

MAX n e = 108 

T M - T i ATP 

°R 

(°K) 

-143.3 

(-79.6) 

-158.7 

(-88.2) 

-1 59.0 
(-88.4) 

-159.0 

(-88.4) 

-159.1 

(-88.5) 

-159.1 

(-88.5) 

-159.1 

(-88.5) 

-159.0 

(-88.4) 


SLAB, 
r M 


Tu = 350°F 


(450 U K) 


I 

TEMPERATURE 
AT P AFTER 
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Baseline Method and Its Computer Time Breakdown 

In order to direct attention to the most time-consuming steps of the calculations, a 
baseline method is first selected and timed to serve as a reference. The baseline meth- 
od consists of standard procedures for the items a) to d) above, namely: 

a) all the E&E’s of each A matrix are calculated 

b) all terms eV are used in all matrix operations 

c) the complete set of E&E’s are obtained via a typical modern transformation 
method for real symmetric matrices, i. e. , an n-step reduction to tri-diagonal 
form followed by the convergent Q-R iteration (ref. 15). Standard routines are 
used (ref. 16). 

d) The analytical expression, eq. (13) is used for the Gjj. 

Table II gives for the baseline method the breakdown of computer time for each of 
the four key steps which constitute one h iteration. This table reflects the medium- 
size problems, n~ 100. The variation with n in the 50 to 100 range is shown in fig. 8. 
Note that in the baseline method, the machine time varies like n 3 and therefore little 
hope can be held that 3D problems can be handled in this way. 

Time Savings Using Dominant E&E’s in the Matrix Operations 

As a first step in the computer-time minimization task, let us evaluate the poten- 
tial of item b). 

It follows from Table II and fig. 8, that no time savings can be accrued in the cal- 
culation of the temperature errors since the computer time is already negligible. There 
remains the possibility of dropping the subdominant-eigenvalues terms in the calcula- 
tion (through formula (13)) of the G { j which takes up most of the time. Unfortunately, 
it turns out that while it is quite possible in each h iterate to drop the subdominant 
terms in the temperature calculation, it is not possible to do so in the Gy calculation 
through formula (13). To see this, let us consider first the temperature calculation in 
each h iterate. 

We have already verified (Table I) that accurate values of the temperature can be 
obtained with only a handful of dominant E&E’s. To examine the problem a little more 
formally, let us derive the same result formally from the equations as follows. 

In each h iteration, for a set of h iterates, the temperatures at the time of melt 
are calculated from 

/ ^ 1,n e- x t 

Ti(t,)=E P u eVi 

i 
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TABLE II. COMPUTER TIME BREAKDOWN IN THE BASELINE NUMERICAL METHOD 


METHOD 

IBM 370/165 (FORTRAN G) TIME FOR ONE ITERATION IN h, SEC. 

STEP 1 

COMPUTE 

E&E'S 

STEP 2 
COMPUTE 
TEMPERATURE 
ERRORS, T m -Tj 

STEP 3 
COMPUTE 
INFLUENCE 
COEFF., Gjj 

STEP 4 

COMPUTE 

UPDATED 

h i 

TOTAL 

BASELINE 
(ALL E&E'S + 
ANALYTIC Gjj) 

-24 

-0 

-79 

-0 

~103 


NOTES: IBM PILOT CODE 

SLAB, HEATED ON 2 SIDES 

TOTAL NUMBER OF ELEMENTS = 108, SURFACE ELEM = 2 X 12 - 24 
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FIGURE 8. MACHINE TIME REQUIRED BY BASELINE VERSION 
OF THE NUMERICAL METHOD 
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where T^) is the temperature at the i-th surface element at the i-th melt time, t t . P 
is a rectangular matrix and the X’s are the eigenvalues, assumed to be ordered so that 
0 \ ^ X 2 ^ A n • Both Pjj and X } depend on h, the heat transfer coefficient for the 

k-th surface element. For any case which involves a nontrivial transient, < 0, and h k > 0 
for at least one k; for the trivial case, h k = 0, all k, Ai=0, Pn-T^O), P^^O, j>l. For 
the latter trivial case clearly there is one dominant eigenvalue, Ai = 0; the practical 
cases of interest involve such relatively small positive values of h k that there is a 
“cluster” of small dominant negative eigenvalues with a complementary cluster of fairly 
large negative eigenvalues. Thus, it usually suffices to sum over the first n e «n terms. 
We shall now derive the explicit dependence of the matrix P on the eigenvectors and 
show that the first n^ columns of P depend only on the n e dominant eigenvectors. The 
equations which T t (t) satisfy are, in partitioned matrix form: 


r— 

— ] 


[— 



— 


— — 


— “1 


— 


0 


T 

= 

A 11 

Aji 


T 

i 

T(0) 

= T 0 

1 

0 

m 2 


t 


A 21 

A 22 


T 


T 


1 


where M*, M 2 are positive diagonal matrices, T is the vector of surface temperatures, 
t that of internal temperature, A^, A 22 symmetric matrices, A^ the transpose of A 21 
and 1 a vector whose components are all ones of appropriate order. The k-th diagonal 
element of the A tl matrix is a linear function of h k . The A matrix is singular for h k = 0 
with 1 an eigenvector. The differential equation can be converted to the form: 
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The solution for T(t) can then be written as: 
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where is a diagonal matrix of dominant eigenvalues of { }, A 2 that of subdominant 
eigenvalues, and the corresponding matrix of eigenvectors is: 


rvn v 12 -i 
Lv 21 v 22 J 

Assuming e A2t is small enough to ignore for t = one of the melt times and considering 
the effect of the matrix multiplications: 

T(t) « [Mj 1 /2 Vn] [e A i t ][v7i Mi 1 /2 1 + vJi M 2 /2 T] = T(t) 
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Now, noticing that the premultiplication of a vector v by a diagonal matrix I D] is equiva - 
lent to the premultiplications of the vector col{[D]} by the diagonal matrix, [diag{v}] the 
result is: 

T (t) = [M } 1 / 2 V U ] [diag {V?i m [ 72 T+ V?i mJ /2 !}] col {[ e^*]} 

= [Pi] col{[e Alt ]} 

where Pi is a function of only [V11/V21I, the dominant eigenvectors. Therefore, we have 
shown even formally that the temperature equation admits an approximate (but as we 
have seen accurate) expression in terms of only the dominant eigenvalues and eigenvec- 
tors.* 


Now the same method used for the temperature errors permits to answer in the 
negative the question as to whether the dominant eigenvalues approximation is accu- 
rate for the calculation of 9 Ti/ 9 h k . 


Consider the following development. 


®TiM_ v 
9h k "t 



+ tiP tj 



e^i 



By similar reasoning the sum can usually be truncated at n e «n terms but the diffi- 
culty is in accurately evaluating R ljk , j $ n e , as a function of only dominant eigen- 
values and eigenvectors. Since R is a third order tensor, the analysis is more in- 
volved . In order to make the point we shall consider the explicit dependence of the 
infinitesimal influence coefficient on the eigenvalues and eigenvectors for the simple 
case: 

i) 2 equal elements, 1 surface element T t and 1 'interior element T 2 , 
with T awi = T a „ 2 = 2 T 0 and uniform initial temperature equal to T 0 

ii) M 1 = M 2 = i 

iii) 1 dominant eigenvalue and 2 subdominant I X t I « I Xj I ; the dominant vector 
is [V11/V21] and the subdominant [V12/V22]. 

The explicit solution of this problem is (temperatures are made nondimensional with 
the initial temperature, times with Ax/«, the heat transfer coefficient with 2 k/ Ax; 

Ax is the size of the element): 


*By implication, then, also the direct heat conduction problem, when solved by the 
eigenvalue method , admits a simplification in terms of only the dominant eigenvalues 
and eigenvectors. 
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It follows: 

r T ‘ ( " ) i=r 2 i . 

Lt 2 ( oo )J I2J 

Let 

hi rv B v 21 -|i-t 1 -t 1 («>)i 
Ur Lv 12 v 22 J Lt 2 — t 2 (°°) J ’ 

where V is orthogonal so 


" T i (t) ~ 

_pi(«n 

J V11 

Vl2 _ 

”Tj(t) ” 

- T 2 (t). 

Lt 2 («)J 

<r 

to 

H-*- 

V 22 - 

-T 2 (t). 


Then 

r*i-|_r*i °ir Ti i r Ti(o) i_r vii+V21 ] 

U 2 J “Lo x 2 JLT 2 J i LT 2 (o)J“ Lv 12 + v 22 J 

Since V can be chosen so that 

rvn v 21 -| r-i-2h 1 -jrv u v«i [A, on 
lv 12 v 22 JL 1 -1 JLv 21 v 22 _tLo xj 


where 

Xj = — (1+h-r); X 2 = — (1 + h + r) ; r = V'l + h 
V 11 =V0.5/[l+h z +hr]; V 21 = V 11 (h + r) 

V 12 =V0.5/U+h z -hr]; V 22 = V 12 (h-r) 

For every h> 0, there is at >0 such that: 
e x iy e *2‘ = e <M-*2)t»i 

for t>t. For those values of h and t, we say that Xj dominates X 2 . It follows that 
T, (t) = 2 — Pj e Xlt -P 2 e Xz * 
where 

Pl = V„(V 11+ V 21 ); P 2 = V 12 (V 12 +V 22 ) 
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Note that the dominant term can be evaluated from knowledge of the dominant eigenval- 
ue, X t , and eigenvector, [Vn/V 2 i], alone. Call this truncated approximation, tj (t) 
where T 1 (t)=2 — P t e Xlt . It is this kind of approximation, retaining perhaps five terms 
out of one hundred, which is used in the difference quotient estimate of the derivative 
of the temperature at a surface element with respect to a heat transfer coefficient 

aT At Tjfh + Ah) — T t (h) AT i Tj (h + Ah) - T t (h) 

9h “ Ah Ah “Ah “ Ah 

Note that there are two levels of error: (1) for the accuracy of the true difference quo- 
tient in estimating the derivative and (2) for the accuracy in the truncated approxima- 
tion to the true difference quotient. For the success of the iteration for estimating heat 
transfer coefficients from melt times, the latter error is surely more serious. This 
is so since generalized secant methods tend to give convergence rates that are compar- 
able to Newton’s method. 


In order to indicate the altogether different nature of the series for the analytical 
formula for dTi/dh, consider the following expression that applies to Eq. (13): 



(t) = (R tl t+ R 12 ) e x l* + (R 22 t - R 12 ) e*2* 


where 

R 1 l = 2Vf 1 <V U + V 2I ) 

Ri2= 2V 11 V 12 (2ViiVi2+ Vi 2 V 2 i + V 11 V 22 )/(X 1 -X2> 

R 22 = 2V 22 (Vi 2 + V 22 ) 

As in the expression for T t (t), the first term surely dominates; however, the coefficient 
has a component, R 12 , which depends on all eigenvalues and eigenvectors. That this is 
generally true follows from a result in matrix perturbation theory (Ref. 20) that the de- 
rivative of any eigenvector depends on all eigenvalues and eigenvectors. R 12 is in fact 
the derivative of P t which itself depends only on the dominant eigenvector, [Vn/V 21 ], 
Define the dominant term in 9T/9h as 9't/9h which is the limit of the dominant term in 
the difference quotient as Ah— -0 
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Note that 6T/6h is not the limit of the truncated difference quotient as Ah-* 0. 

In order to illustrate the effectiveness of the truncated difference quotient by com- 
parison with the truncated analytical derivatives, expand the above expressions in terms 
of h and t. 


t)= 2_r°-_5 ( v h+r) i 

’ 1 2 [ l + h 2 + hr J 

t) = 2 -[ M l V h+r) 1 

’ * |_ l + h^ + hr J 


e -(Uh-r)t 

e -(l+b-r)t 


0.5(1 +h — r)~ l _( 1+hfr) t 

1 + h 2 — hr J e 


AT T(h + Ah)-T(h) 
Ah Ah 


6T 

6h 

0T 

8h 

9T 

9h 


•(l*h-r)t 


t) (h, t) + 0. 5 e ‘ 


(l+h-r)t 



e -(l + hfr)t 


where 

r= Vl+h 2 . 

a 

Figure 9 illustrates that the truncated difference quotient, AT/ Ah, is a very good 
estimate of the true derivative, 9 T/9h, for values of h and t which satisfy the domi- 
nance condition that T is a very good approximation to T. The data for the illustration 
were chosen to correspond to a typical phase change test condition, i.e., T lnlt = 76°F 
(298 °K), T aw = 946°F (781 °K), T me u = 300 °F (422 °K), Wt=2sec. The value used for 
Ah ensures that AT/Ah estimates 9T/9h within 0.5%, while I AT I remains larger than 
0.1%. The latter condition guarantees that the computed version of A^/Ah will have 
three correct digits, for example, if the computed version of T has six correct digits. 

The conclusion then is that item b), dropping the subdominant eigenvalues terms 
in the calculation of the temperature errors and the influence coefficient Gjj, does not 
offer potential for computer time savings since the temperature calculation takes up 
negligible time and the G t j expression cannot be cast approximately in terms of only 
dominant eigenvalues. 


Time Savings by Choosing an Appropriate Method for Calculation of the Eigenvalues 

and Eigenvectors 


Next let us consider item c). The baseline method as indicated above uses, for 
the E&E’s, a transformation-type method for real symmetric matrices. But an attrac- 
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tive alternative should be an iteration-type method because good starting approxima- 
tion for the E&E’s are available from the previous h iterations. Naturally a method to 
obtain a starting approximation must be provided for the zeroth h iteration, but this is 
less crucial a step since it occurs only once and the really important time savings 
should come from all the subsequent h iterations. The best iteration-type method for 
our application seems to be the Jennings algorithm /Ref. 17 & 18) since (i) it can take 
advantage easily of the sparseness and bandedness of the A matrix, (ii) it is just about 
the only method for dominant E&E’s, (iii) it is very efficient if the starting approxima- 
tion is good — which is a very useful feature in our case when the h iterations provide 
us with increasingly better starting approximations of the E&E’s. 

Therefore, a specially adapted version of Jennings’ algorithm has been used . Due 
to the distribution of eigenvalues for the heat equation, the implicit form of the inverse 
matrix is used for the eigenvector recursion or power step. The Jennings’ algorithm 
used does not include a refinement, suggested by Clint and Jennings (ref. 18), called 
“Jacobi eigenvalue reduction. ” 

Jennings’ method has been adapted by coding the matrix multiplications to take full 
advantage of the sparse and banded form of the real symmetric matrix which arises in 
this problem. Multiplications which would produce zero results are thus omitted. This 
should result in a considerable time savings in the eigenvalue-eigenvector refinement 
over that of the standard Jennings’ algorithm as coded by Vachris (ref. 19), for ex- 
ample. The Jennings’ algorithm turns out to produce considerable time savings com- 
pared to the base line method as fig. 10 shows. 

The semi-discrete form of the heat equation results in a matrix for which the most 
important eigenvalues are the smallest ones. Jennings’ method requires dominant 
eigenvalues to be the largest ones. There are two transformations which could ac- 
complish this. The simpler one to apply involves shifting all the eigenvalues by a uni- 
form amount. This, unfortunately, results in inherently slow convergence of Jennings’ 
method since the ratio of the largest unimportant (shifted) eigenvalue to the smallest 
important one is around 0.99. The approach which we have implemented results in 
faster convergence with values of the above ratio of around 0.4. This method utilizes 
a special efficient algorithm for inverting a banded, positive definite, symmetric ma- 
trix and is fully described in ref. 20. 

In conclusion, item c, the use of an optimum method such as the Jennings’ algorithm 
provides interesting time savings when only the dominant E&E’s are required. Simul- 
taneously, it follows also that item a), whether it is possible to save machine time by 
not calculating the subdominant eigenvalues, is answered in the affirmative. 
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FIGURE 10. MACHINE TIME REQUIRED BY TWO OF THE BEST METHODS FOR ALL 
OR THE DOMINANT EIGENVALUES AND EIGENVECTORS 
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Incidentally, a different adaptation of Jennings’ algorithm to problems similar to 
our is reported by Rutishauser (ref. 21 and 22).* However, Rutishauser’s adaptation 
is very complicated and contains various sophisticated features. The machine time 
quoted for a 70 element problem should correspond to some 30 sec. on the IBM 370/165. 
Naturally, if this adaptation were applied to our problem, it would result in running 
times that are larger than the QR algorithm in the baseline - method . 

Minimum Computer Time Alternatives 

Summarizing the discussion up to now, we have eliminated item b) and we have se- 
lected the optimum method for item c). Therefore, the alternatives remaining are a), 
whether to calculate all or dominant eigenvalues, and d) whether to calculate G 13 analyt- 
ically or numerically. These alternatives form the matrix given in Table III. The main 
conclusion so far is that since the Gy take up most of the computer time in the baseline 
method, minimum computer time will be achieved by (i) calculating only the dominant 
eigenvalues; (ii) calculating the G tJ through the numerical method. Indeed, timing of the 
alternatives gives the results already indicated in Table III. 

The machine time required by the numerical G tj and the dominant eigenvalues de- 
pends upon the number of dominant eigenvalues and naturally decreases drastically with 
that number. The optimum number is the one that minimizes the CPU time maintaining 
a good accuracy in the results. As expected, it turns out that this number is very 
small (fig. 11) and the number of iterations unchanged at constant accuracy (fig. 12). 

It is veiy satisfying that both accuracy and the number of iteration remains substan- 
tially constant for a wide range of attractive values, so that a universal working value 
can be easily selected. 


♦Initially, Bauer (ref. 23) generalized Jennings basic idea to nonsymmetric matrices. 
Rutishauser (ref. 21 & 22) subsequently specialized and refined what he refers to as 
Bauer’s simultaneous iteration method to the symmetric case. Rutishauser (ref. 22, 
p. 221) discusses an example which has all of the properties of our problem. The 
system matrix is banded, sparse, and definite, the dominant eigenvalues are those 
which are smaller in absolute value and the relative magnitudes of the shifted eigenval- 
ues to their now smaller immediate neighbors is very close to 1. For Rutishauser’s 
example the ratio is 0.999 while for our example it is 0.99. The shifting is required 
in order that the working matrix have, as its dominant eigenvalues, the largest ones 
in absolute value. Of course, the alternative to shifting is inverting by the Cholesky 
factorization method. The expense of the factorization and subsequent multiplications 
by now dense, banded matrices in Rutishauser’s case more than compensated for the 
slow convergence and poor final accuracy of the shifted result. 
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TABLE III. VARIANTS OF THE NUMERICAL METHOD THAT HAVE BEEN STUDIED 
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FIGURE 12. NUMBER OF ITERATIONS AND ACCURACY USING A LOW NUMBER OF DOMINANT 
EIGENVALUES AND NUMERICAL INFLUENCE COEFFICIENTS g n 





The machine time required by the alternative A2 a in the matrix (Table III) de- 
pends also on the number of dominant E&E’s selected. Again the selection can be 
optimized, see fig. 13 and 14. 

The machine time required by the alternative A2/3) is the same as B2/3) since the 
difference — for the same number of dominant eigenvalues — is just whether or not all 
the E&E’s are used in the matrix operations and we have already seen that the differ- 
ence is negligible. 

The optimum machine time in each alternative is compared in fig. 15. The con- 
clusion is that numerical and only dominant E&E’s represent the minimum machine 
time alternative and therefore the one used for the operational computer code. The 
variation with the number of elements of the optimum alternative is shown in fig. 16 . 
The important result in this figure is that the selected alternative not only has the low- 
est machine time, but also its machine time varies as n 2 against approximately n 3 of 
the baseline method. This is a remarkable feature of the method developed here that 
permits to run relatively large problems. 
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FIGURE 14. NUMBER OF ITERATIONS AND ACCURACY USING A LOW NUMBER OF DOMINANT 
EIGENVALUES AND ANALYTIC INFLUENCE COEFFICIENTS Gj: 
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FIGURE 15. COMPARISON OF THE SMALLEST COMPUTER TIME OF EACH 
VARIANT OF THE NUMERICAL METHOD 







Accuracy and Computer times of the Numerical Method Developed 


The method presented, even when run with a minimal number of dominant E&E’s, 
is very accurate. This is brought out in a typical case by fig. 17, where the results of 
the same problem are presented with increasing number of points. 

Note that the three values used in fig. 17 for the total number of points range from 
a very low number (n ~ 30), to an average number (n ~ 100). The experience to date 
suggest that n ~ 100 is more than adequate for typical cases and target accuracies of hj 
of 1 to 2%. Unfortunately, it is not easy to compare the accuracy of this method with 
others, because no other inverse method (with lateral conduction) appears to be avail- 
able and nor any exact analytical solution has been tabulated. 

A strong indication of the soundness of the method is the manner of convergence 
during the hj iteration. Fig. 18 shows the decrease in error (from a final result) at 
each iteration. 

As far as tolerances used to declare the iteration converged, there is no point in 
requiring very tight tolerances in light of the soundness of the method. Fig. 19 shows 
that with 0.5-1% error, the number of iterations can be reduced to typically 4 for slab- 
like problems and 2 for l.e. special problems. 

The computer code has been checked out by running numerous test cases. In this 
section we show the results for the sample problem supplied by NASA Langley Research 
Center. The computed heat transfer coefficients, fig. 19a, match the original values 
within ± 5%. For this case the phase-change temperature was considered variable and 
the time of melt was a constant 3.6 seconds. The code was also checked in the more 
conventional mode of constant phase-change temperature with a variable time of melt. 

To obtain the necessary input information, the NASA supplied run was first duplicated 
on our thermal analyzer (direct problem). By plotting the temperature histories, the 
times of melt corresponding to T cp =1000°R were read off this curve. This information, 
along with the resulting h values from the code is shown in figure 19b. The accuracy 
in h is again within ± 5% . 

A quantitative indication of the potential machine time reduction exploited during 
the optimization effort carried out in this study is given by fig. 20. This is a measure 
of the method’s potential that was announced at the outset. 

Finally, a point of the maximum importance is the fact that this method appears to 
require machine time proportional to n 2 (n is the number of elements) rather than n 3 
that appears unavoidable for implicit methods with time discretized. This is very im- 
portant if one insists on or needs to handle problems in the n ~ 200 range or above. 
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FIGURE 17. ACCURACY OF HEAT TRANSFER COEFFICIENT VERSUS 
NUMBER OF ELtMENTS 
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FIGURE 18. TYPICAL BEHAVIOR OF THE HEAT TRANSFER COEFFICIENT 
DURING CONVERGENCE TO SOLUTION 
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FIGURE 19. TYPICAL NUMBER OF ITERATIONS REQUIRED BY A 
GIVEN TOLERANCE IN THE ITERATION 
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Prospectives for Further Improvements 


In spite of the considerable improvements realized during the optimization effort 
carried out, it is quite likely that considerable more savings on machine time and stor- 
age can be obtained. 

First of all, the Gjj calculation seems not to be necessary for every hj iteration. 
Just doing it every two hj iteration, will cut the machine time by some 40% (assuming 
even number of iterations and the machine time breakdown per iteration shown in table 

n). 

Moreover, the Newton iteration is the simplest and most straightforward procedure, 
but it is probable that a more sophisticated procedure may cut the machine time per 
iteration. This should be worth another 25% savings. 

Finally, the present method could be applied also to direct heat transfer problems. 
Of course, there would be a considerable reduction in computer time over the inverse 
method. A rough estimate of this reduction is as follows. Currently, a typical slab- 
like inverse problem uses about (k + 1) I calculations of E&E’s where k is the number 
of surface elements and I is the number of iterations. A direct problem requires just 
one calculation of E&E’s. However, the calculations of E&E’s in inverse problems ac- 
celerate as convergence is approached; this is particularly true of the k calculations 
needed for the Gjj . Probably then the ratio CPU time, direct versus inverse problem, 
is more likely equal to the ratio of iterations in E&E’s calculations via the Jennings 
algorithm or i 2 / (ki 0 + i 2 ) I where i 0 is the number of Jennings iterations during the G tj 
calculation, ij that for each I iteration and i 2 that for a calculation of the E&E’s starting 
with poor zeroth-order guesses. Typical values found have been: i 0 = 2, q = 4, i 2 = 20, 

1 = 4. Therefore, the expectation is that for a typical slab-like 2-sided problem where 
k ~ 20, the direct method will require — with the method as it stands — only some 11% 
of the CPU time of the corresponding inverse problem. 

TYPICAL RESULTS ON LATERAL CONDUCTION EFFECTS 

The computer tool developed (the code is described in the appendices; its name is 
CAPE (Conduction Analysis Program using Eigenvalues) ) makes it possible to obtain 
quantitative results on lateral conduction effects on the data reduction. While the em- 
phasis in this study was on developing the computer tool, in this section we briefly pre- 
sent some of the typical results that have been obtained. 

The first question is naturally, how large are the differences in the h obtained with 
lateral conduction and without? Figure 21 gives the comparison in a representative 
case, a slab of the small size found for example on the fin of orbiter models of 1 ft. 
length. As one would expect, lateral conduction returns h’s with higher peaks. 
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FIGURE 21. EFFECT OF LATERAL CONDUCTION IN THE DATA REDUCTION 
FOR A SLAB OF INFINITE DEPTH 
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FIGURE 22. EXAMPLE OF COMPLEX INTERPLAY BETWEEN LATERAL CONDUCTION AND FINITE 
THICKNESS EFFECTS IN SLAB HEATED FROM TWO SIDES 


Lateral conduction plays a somewhat more complicated role in the vicinity of a 
large h(x) gradient, when the slab is thin and in fact heated on both sides as in figure 
22 . 

When the h(x) gradient is not so large, the effects are smaller and qualitatively ob- 
vious, see figure 23 where the finite thickness is obtained via the charts of Hunt et al 
(ref. 4) 

In wing l.e. ’s of the size typical of orbiter models, the lateral conduction and finite- 
thickness effects are both large. It is no surprise that the results are rather different 
than by semi-infinite slab reduction, the only procedure available before this study. 
Figure 24 gives a good idea of the typical errors one would incur if one were to neglect 
lateral conduction and finite thickness in such small models. * In this case, it is not 
obvious how to separate the two effects, lateral conduction from finite thickness, in the 
nose region. We believe that the basic reason for the large differences in h over the 
wedge portion is the finite thickness effect. This effect on the wedge portion is very 
important as at/1 2 is roughly 3. 5 and the yes-no chart indicates that corrections are 
definitely needed. A finite slab calculation over the wedge should eliminate a major 
portion, but not the entire discrepancy. Lateral conduction effects should still be non- 
negligible over the wedge. An assessment is obtained by considering that the heat input 
near the nose reaches the wedge within the times of melt, as the length reached is about 
0.017 ft and the wedge portion extends from 0.004 to 0.008 ft. 

Radiation, as expected, is negligible in typical situations involving Stycast. Fig. 

25 gives typical quantitative results on this matter. 

Properties variation with temperature in this case of Stycast are also negligible, 
see fig. 26, which represents ar. extreme, if somewhat artificial case, in that the melt 
temperatures were relatively high and the c p variation with temperature was somewhat 
pronounced, simulating a behavior found in some previously heated samples of pres- 
sure-cast Stycast 2762. The properties used were measured by Revenko and Hansen 
of the Grumman Aerospace Corporation. 

Finally, in exercising this data-reduction computer tool in the presence of signifi- 
cant lateral conduction, we found sometimes large sensitivity of the h’s from the (in- 
evitable) errors in measuring the melt times. This problem, that turns out to be a ba- 
sic problem of the phase-change-paint technique in the presence of either finite thick- 
ness or lateral conduction, is discussed separately below. 


♦Incidentally, it should be kept in mind that such lateral conduction errors are not pecu- 
liar to the phase-change technique, as severe effects are also found with thin-skin ther- 
mocouples in regions of large lateral heating gradients. 
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Fig. 23 Example of Simple Interplay Between Lateral Conduction and Finite Thickness 
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WHEN IS THE COMPUTER TOOL NECESSARY 
TO REDUCE THE DATA? 

Parameters Characterizing the 
Correction to ‘Semi-Infinite-Slab ’ Formula 

The problem, of course, is one of providing a single set of “yes-no” charts that 
will quickly tell the test engineer whether or not he needs to reduce the data via the 
computer tool that accounts for lateral conduction. More generally, the “yes-no” chart 
should tell whether the error over the ‘semi-infinite slab’ reduction (ref. 1) becomes 
excessive, say 10%, because of (1) lateral conduction, (2) surface curvature, and (3) 
finite thickness. 

In order to maintain simplicity and generality in the chart, we isolated four non- 
dimensional parameters that globally characterize the three effects. The approach is 
as follows. 

The criterion on whether a slab is thick or thin is given by the approximate rule 
(Jones and Hunt, ref. 1) . 

a t/\2 _ H 0.2 semi-infinite slab 
1>0.2 thin slab 

If the slab is thin the corrections of Hunt and Pitts (ref. 3,4) are required. This crite- 
rion, of course, is strictly valid only when other corrections are negligible, but it in- 
dicates the key parameter characterizing thickness corrections. 

A criterion for curvature can be readily obtained by comparing transient heating of 
slabs and cylinders under otherwise comparable conditions. The solutions for these 
geometries are available in the standard heat transfer literature for many types of 
boundary conditions, (e.g. , ref. 24, fig. 4-6 and 4-7). Using the solutions for given 
and constant heat flux, it turns out that the difference in heat transfer coefficients for 
the slab and cylinder exceeds 10% when at/R 2 c is greater than 0.12. This establishes 
the limit 

, , (< 0.12 no curvature correction 

at/Rl ^ 

(>0.12 curvature corrections needed 

This identifies a synthetic non-dimensional parameter for curvature corrections. The 
parameter has the obvious meaning of the distance of heat propagation compared to rele- 
vant length, the radius of curvature. Incidentally, at first glance it might appear that 
curvature corrections are more restrictive than thickness corrections. Actually, this 
is rarely the case. Except for the leading edges, the curvature is generally much 
greater than the thickness. Consequently, oit/R\ would have a low value and curvature 
corrections would be negligible. 
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The situation of the lateral conduction is a little more complicated. The gradients 
in the heat-transfer coefficients can be essentially of two types, as shown in the sketch 
of fig. 27. When the gradients are very steep, we have a ‘step’ or a ‘top-hat’. These 
two cases, even though severe, represent a useful reference. 

In seeking synthetic parameters, analytical solutions to the problem of variable 
heat flux again provide a precious guide. Analytical solutions for a ‘top-hat’ of given 
heat flux, are available for both the two-dimensional (strip) and three dimensional 
(spot) cases on a semi- infinite slab (ref. 7). The corresponding temperature distribu- 
tions are reinterpreted in terms of required corrections to the semi-infinite slab heat 
transfer coefficient in fig. 28. It is seen that the effect of lateral conduction ofi the 
heat transfer coefficient is “correlatable” by the single parameter a// at. This param- 
eter can be put in the form at/a 2 completely similar to the parameter characterizing 
the thickness of the slab. The only difference is that the length scale is the width of the 
“top-hat” instead of the thickness of the slab. As might be expected, the lateral conduc- 
tion effects for the spot are greater than for the strip*. If one takes errors in q at the 
center of the hat in excess of 10% as significant, one can determine the values of the 
a/ Vert parameter at which this occurs: 

a/ Vat at/a 2 

strip 1.5 0.445 

spot 2.0 0.25 

The parameter at/a 2 characterizes the errors, since the point at which one is most 
interested to find h is precisely the center of the hat. In other words, the distribution 
of the heat input carries explicitly the relevant length scale. 

The situation is a little different for the step distribution because here there is no 
length associated with the heat input. This is a reminder that in the immediate region 
of the discontinuity, lateral conduction is always important. But as the discontinuity is 
a gross schematization of a large gradient, the existence pf lateral conduction at the 
discontinuity is just a defect of the schematization. Appropriately, in the case of the 
‘top-hat’, the position for judging whether or not there are significant lateral conduction 
effects was not at the discontinuity, but at the center of the hat under the presumption 
that the length a is not much smaller than the length scales over which the problem is 
examined. It follows that the meaning of the parameter emerging (from fig. 28) is the 
usual one, a diffusion length Vot/l.45 becoming about equal to the semi -width of the 
hat and thereby affecting the conditions at the center of the hat. In the same fashion, 
for the step, lateral conduction will be important at distances less than Vot/l.45. The 
question is whether such distances are of the order we are considering. If the slab 
is of width w, lengths about O.lw must surely be of interest. In other words, some 
ten time-of-melt measurements are conceivable over the width w. This fixes the 


*The spot problem is a xi symmetric and cannot be handled by the code as it stands. 
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FIGURE 27. SKETCH OF THE TWO ESSENTIAL TYPES OF h DISTRIBUTION 





minimum length scale of interest. Therefore, to be noticeable, lateral conduction 
in the case of the hat must diffuse at least 


Vat 

1.45 


)0,lw 


or 

10 2 ^ 0.47 

This is the appropriate parameter for step distributions. 


‘Yes-No’ Charts 

The limits for lateral conduction can now be combined with the criteria for thick- 
ness and curvature in a single “yes-no” chart, fig. 29. This chart maps out the region 
where different corrections are needed. As special case the chart incorporates the 
criterion for thickness in the absence of lateral conduction and curvature (Jones and 
Hunt, ref. 1) and therefore indicates the ranges where the corrections of Hunt and Pitts 
(ref. 4) maybe applied. 

The yes-no chart should also be used in deciding where to place the adiabatic 
boundary in l.e. problems. The distance of this boundary from the stagnation point 
should at least be greater than the value obtained by at/1 2 = 0. 2. The results of fig. 24 
do not satisfy this criterion by a factor of 2. 
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ACCURACY OF THE PHASE-CHANGE 
PAINT TECHNIQUE IN THE 
PRESENCE OF FINITE THICKNESS 
AND LATERAL CONDUCTION EFFECTS 

Effects of Inaccuracies in the Tunnel Input 

In a comparison between results from different computer tools such as fig. 30, 
where the direct and inverse problems were run with two different computer codes, 
some waviness was observed in the h distribution from CAPE. We have never encoun- 
tered this behavior when CAPE was used to solve both the direct and inverse problems. 
One can speculate that the waviness results from slight inaccuracies in plotting the tem- 
perature, cross plotting the times of melt and reading off the appropriate values for in- 
put to CAPE. Since these functions were performed by eye, some loss in accuracy is 
inevitable. If we deliberately perturb the time of melt distribution, in a simulated data 
reduction, we find the same behavior, as indicated in fig. 31. These perturbations are 
imperceptible on the time of melt curve; however, they resulted in pushing the h distribu- 
tion outside the ±5% error limits. Thus, any waviness in the time of melt distribution 
appears to be greatly amplified in the h distribution. This is particularly serious where 
points are closely spaced. The same behavior is observed if, rather than controlled 
perturbations of the melt times at specific locations, random melt time perturbations 
of given rms are generated by computer as in fig. 32. 

An essential point that is to be stressed is that this happens only in regions where 
lateral conduction effects are important or, more precisely, when dh/dx are large 
enough to make them so. See, in fig. 32, how well h is recovered when dh/dx is not 
inordinately large. The problem disappears, if the melt times are reduced with the 
semi-infinite slab approximation - except of course that the h estimates are wrong. The 
‘waving’ of the solution, and particularly the higher amplitudes obtained by refining the 
grid have a simple physical interpretation. 

One can readily imagine what would occur when the time of melt distribution is ob- 
tained from tunnel measurements with all the inherent uncertainties. Therefore, the 
problem cannot be ignored. In fact, it appears of the greatest importance if one is to 
be able to reduce the tunnel data. 

Smoothing of Tunnel Inputs 

The natural solution to the problem is to 'smooth' the raw time-of-melt data prior 
to use in CAPE or any other inverse tool, at least when lateral conduction is non-negli- 
gible and is to be taken into account. Naturally, smoothing of raw-data is nothing new 
in wind-tunnel data reduction and is included in many automatic data reduction systems. 

As tool for smoothing ‘noisy’ melt time data, we can borrow a recently developed 
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DISTANCE 


FIGURE 30. WIGGLES GENERATED IN THE INVERSE PROBLEM WHEN THE INPUT DATA 
ARE OBTAINED BY ANOTHER HEAT CONDUCTION CODE 
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h, BTU/sec ft z °F TIME OF MELT, SEC 



FIGURE 31. SENSITIVITY OF INVERSE SOLUTION TO DELIBERATE 
PERTURBATION OF THE TIMES OF MELT 




FIG. 32 - DATA REDUCTION WITH RANDOMLY PERTURBED MELT TIMES 
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h, W/m 2 


technique, the ‘smoothing spline* (ref. 25, 26). In essence, the method allows one to 
choose an arbitrary degree of smoothness by specifying a parameter R. When R = 0 
there is no smoothing and the method produces a spline fit through all the data. When 
R = 1 one obtains ultimate smoothing, i. e. , a straight line. For intermediate values of 
R, fig. 33 shows some examples of various degrees of smoothing. The right amount of 
smoothing is to be found by trial and error. 

In detail, the smoothing spline consists of the following. Given a set of raw data 
g t associated with a coordinate s lf i = m to n, in the interval for s{0, l}, one computes 
a set of f t (s), a third-order spline for each interval. The four parameters in each f^s) 
are chosen so as to minimize the expression: 



where R t = 1 - R 2 , R 2 = (1 — R) 4 and Q are chosen parameters that control the degree of 
smoothing. R 2 = 0 gives the least square straight line fit, while Rj = 0 gives the normal 
spline. 

The technique has been applied to the randomly perturbed times of melt of fig. 32. 
The results were not favorable as shown in fig. 34. Indeed the smoothing process re- 
duces the waviness in the h results, but this does not prevent unacceptable deviations 
of the h values near the peak. The degree of smoothness infig. 34, R^0.3, is justabout 
the maximum values that one can reasonably use, since as R increases the times of 
melt are modified to the point of changing the problem, see fig. 35. Naturally, in going 
from R = 0.1 to R = 0.3 there are definite improvements, as fig. 36 shows, but the re- 
sults remain unsatisfactory. 

If, instead of smoothing, the melt times are just fitted and the fit is used in the data 
reduction, the results are equally unsatisfying, apparently because the fit is poor at the 
ends of the interval. This is a facet of the problem that should be looked at more close- 
ly as perhaps there is a chance for at least better results than fig. 36. 

It becomes clear that the problem is of deeper nature and is not solved by straight- 
forward treatment of the raw data before reduction. 

Interpretation of the Difficulties Encountered 

The difficulties experienced lead us to explore the fundamental problem of the sen- 
sitivity of the phase-change paint technique to errors in the times of melt. Clearly, this 
is the crucial point emerged above, i. e. , the fact that small variations in the times of 
melt give large variation of h resulting from the inverse problem. 

A basic condition for a sound experimental technique requires that the errors in the 
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SMOOTHING SPLINE 




FIGURE 33. SMOOTHING SPLINE FOR TREATMENT OF "NOISY" EXPERIMENTAL DATA 
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FIGURE 34. DATA REDUCTION AFTER STRAIGHTFORWARD SMOOTHING OF THE MELT TIMES 
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FIGURE 35. EXCESSIVE SMOOTHING MODIFIES INTRINSIC CHARACTERISTICS 
OF GIVEN SET OF MELT TIMES. 
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measured t,,, are not magnified into much larger errors of the results of h, i.e. , the 
error amplification 


T 5T =0(1) 

h t 


so that the % error in hj is about the % error in the t t . As a typical value for accurate 
and sound techniques, we can take error amplification less or equal to 2. 


The question is what values have the (tj/h t ) (dhj/Btj) in the typical applications? The 
sensitivity found suggests that the (t^/Tij ) (dhi/Btj-) are unduly large. This can be checked 
by calculating them. 


Error Amplification in the 
Semi-Infinite Slab Reduction 

In situations where the semi-infinite slab reduction is accurate, the calculation of 
the error amplification is straightforward. By differentiating the explicit expression 
for h (t^ (ref. 9) it follows that 

_ tm,j 3hi -5 for i = j 
h t 9^, 0 for i * j 

Therefore, as one would have suspected from the general success of the phase-change 
paint technique, the technique is sound and the errors are not amplified and most ob- 
viously, are not propagated. 

Error Amplification in the 
Presence of Finite-Thickness Effects 

In the presence of important finite thickness effects, but no lateral conduction, it 
is obvious that 

tnU 9h l 
hi 

if j *i and not on the ‘other side’ of i. Therefore, at each surface point i, there are 
two influence coefficients (calling 1 and 2 the two sides of the slab): 


An 


_ ^ml 
hi 


8h, 

9W 


A 9hl 

12 h, dt^ 
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Naturally, the same considerations are valid for the c other side’ of point i. These co- 
efficients must be obtained numerically, for example with CAPE or with the tool de- 
scribed in ref. 4. This has been done for all slabs heated on two sides with same time 
of melt and same T aw on both sides. In this case there is only one independent non-di- 
mensional parameter, r ml = a t^/l 2 . Therefore, A u and A 12 depend only on this param- 
eter. Of course, as r m i — 0, we recover the semi-infinite slab and therefore A u — -i 
and A 12 “* 0. For arbitrary T m i 5 A n and A 12 are the universal functions given in fig. 37. 

The important result in fig. 37 is that at large r ml A u and A 12 are very large. 
Therefore, the waviness encountered in the data reduction is the result of large error 
amplification intrinsic in the physics of the dependence of the h on the tm when finite- 
thickness effects are important. 

There follows the conclusion that at very large t m and in the absence of lateral con- 
duction, the phase-change paint technique is unsound since it does not satisfy the con- 
dition that A tl and A 12 be less than, say, 2. Moreover, using this limit, one can deduce 
from fig. 37 that the range of applicability in the presence of finite-thickness effects is 
below about T m = 2 . 5 . 

One relatively common situation in which one may encounter r m of the order of 2 is 
the trailing edge of a thin fin section such as the one studied in ref. 4. 

Error Amplification in the 
Presence of Lateral Conduction. 

Strong lateral conduction must also result in large error amplification. This is 
suggested by the Svaviness* experienced in the data reduction process and, more fun- 
damentally, by the physical idea behind the error amplification. Therefore, along lines 
similar to the finite-thickness case, one can determine the limits of applicability of the 
phase-change paint technique in the presence of strong lateral conduction. Also, it 
seems probable that the potential of stainless steel models or data reduction on the pres- 
ence of fully three-dimensional lateral conduction should be seriously questioned. 

These matters have not been looked into, but deserve a close quantitative study. 
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FIGURE 37. AMPLIFICATION FACTORS AND RANGE OF GOOD DATA REDUCTION FOR ONE 
DIMENSIONAL FINITE-THICKNESS SLAB WITH EQUAL TIMES OF MELT 





CONCLUSIONS 


The main conclusions of this study are: 

1. A single and synthetic chart has been developed that permits one to quickly decide 
when lateral conduction and finite-thickness effects, have to be accounted for in the 
data reduction for the three geometries of interest. 

2. A computer tool has been developed to carry out the data reduction for slabs, l.e. 
and arbitrary two-dimensional geometries. For typical slab problems, the ma- 
chine time is about 6 min. on the CDC 6600, with an accuracy of 1 to 2% on the h 
calculation. The labor needed for the preparation of the input data is 15 to 30 min- 
utes per case. 

3. The computer tool can also handle a special data reduction problem for l.e.’s, 
where only one time of melt is supplied. It is now possible to evaluate, by appro- 
priate experimentation, this type of data reduction. The machine time for this 
problem is about 2 min. on the CDC 6600. 

4. Lateral conduction is found to have considerable effects on the h value returned es- 
pecially for l.e. of wings and fins of sizes typical 1 ft models of complete vehicles. 

5. The inrccuracies in time-of-melt data produce considerably amplified inaccuracies 
in .' o a when strong finite-thickness or lateral conduction effects are present. 
Straightforward smoothing of the raw data does not help. This behavior appears to 
limit the range of applicability of the phase change paint technique. For finite- 
thickness effects, the limit is roughly nt m /l 2 ~ 2. 5 where a is the thermal diffusiv- 
ity, t m the time of melt and 1 the thickness of the slab. 

6. As byproduct of the effort, an interesting numerical method has been developed that 
applies to heat conduction ideas evolved in control theory and structural analysis. 
The method can be applied also to direct problems of unsteady heat conduction with 
rather low computer time estimated, in two dimensions and with constant proper- 
ties, at about 11% of the corresponding inverse problem. 
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APPENDIX A 


USER ORIENTED DOCUMENTATION OF THE CODE 

This appendix contains the amount of information strictly necessary to the engi- 
neer to be able to use the code. The details of the code structure are set out in Appen- 
dix B. The code named CAPE (Conduction Analysis Program using Eigenvalues) solves 
the two type of inverse problem mentioned in the introduction, using the method de- 
scribed in the body of the report. CAPE is programmed in Fortran IV for the CDC 6600. 
A version with a few double precision statements is also available for the IBM 360/175. 

The geometries handled by the code are those indicated in Figures 4, 5 and 6. For 
slab-like geometries, CAPE can also handle one-dimensional finite slabs, quasi-one 
dimensional arbitrary geometry cases and semi-infinite (in depth) two-dimensional ge- 
ometries. The calculation of h via the ‘semi -infinite-slab’ (SIS) formulae is also done 
in every case and printed out. 

The grid lay-outs are given in figures 4, 5 and 6. In order to operate the code, 
the first step is to assemble the input information, i.e. , 

a) establish dimensions of the slab, afterbody angle for the l.e. , and geometry 
definition of the arbitrary four-sided geometry; 

b) decide for one of the two types of problems, the one where all the times of 
melt are given or, for the l.e. , the special option where only the minimum 
melt time around the 1. e. is given; establish, for slab-like problems, whether 
both sides are heated or just one. 

c) tabulate the time(s) of melt together with the position on the surface (s) at 
which they are known. 

d) secure the model material properties, i.e. density, conductivity and specific 
heat, whether variable with temperature or not; if the material is Stycast, 
CAPE can be instructed to automatically select constant, but appropriately 
chosen properties. 

e) tabulate the adiabatic wall temperature on the model surface (s), whether con- 
stant or not. 

f) tabulate the melt temperature (s), whether variable on the surface (s) or not. 

g) if a l.e. geometry uses the option of a single time of melt given, obtain the 
tunnel gas conditions, M^,^, 7, p*, T tot , R, Pr, a eff (the tunnel gas is 
treated as having constant: 7, gas constant, Prandtl number, etc.). 

and ot efi are the effective values in cuts normal to the 1. e. 

h) decide whether or not the experimental times of melt need ‘smoothing’ of the 
inevitable scatter. 
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The input cards that need to be prepared for one case are described in table A-l 
together with the definition of the symbols in table A-2 and Fig. A-l, A-2 and A-3. 

The output printed out by CAPE for each case is described in table A-3 (the input 
data printout and the initial calculations), and in table A-4 (the results). This is the 
normal printout containing the information needed by the user. 

Cases can be stacked at will. 

CAPE uses no tapes and no disks. It does use a single overlay consisting of a 
root region and three primary levels. The standard score 3.2 overlay feature is used 
to prepare the overlay. Three additional main programs control each of the primary 
levels, the names being MPCP, MDETRAD and MOUT. CAPE is made up of the main 
program and 33 routines. 

Minimum storage required depends upon the size of the problem, i.e. , the number 
of elements of the problem. This is given in figure A-4. Typically, good accuracy is 
obtained with N ~ 100, for example, N = Lx M = 10 x 9 = 90 for the 1. e. , and N = MxL 
= 12x9 = 108 for a slab-like geometry. For this problem, the storage required by CAPE 
is 60 8 k. CAPE includes a feature that enables the user to keep the storage to a minimum 
for the N selected without having to change dimension cards. Appropriate dimensions 
are set automatically by the code but the total length in decimals has to be input into the 
code. The length can be read from fig. A-4 and the following two cards of the main 
program have to be set as follows: 

DIMENSION S(L) 

CALL SIZE (S, L) 

where L is the length in decimals read from fig. A-4. Naturally, if one is not interested 
in the storage savings that go with this feature, one can just set L, once and for all, for 
the typical N~100 problems. When stacking cases, L must be chosen so as to accom- 
modate the largest case of the stack. 

A guide on the computer time required on the CDC 6600 is given in fig. A-5. Nat- 
urally, the computer time depends somewhat on the problem, so that deviations from 
fig. A-5 are possible. 

Some advice is appropriate on the choice of a few parameters: 

1) for the total number of points, it seems that N~100 gives accuracies of the 
order of 1 or 2% in h and is therefore considered satisfactory since it is about 
one order of magnitude more accurate than the final result, i. e. , the experi- 
mental h. 

2) In order to maintain reasonable balance in accuracy between the two directions 
in a siab, typically one should use numbers of elements LXM such as 9 (in 
depth)xl2 (on the surface), 11x19, etc. Analogously, for the l.e., typically 
MCAP = 10, M = 9, L = 10 (see figures A-l, A-2 and A-3). 
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TABLE A'l. INPUT DATA FORMAT FOR CAPE 


Card No. 

Data 

Format 

Remarks 

1 

NSIDES, L, M, NE, LLL 

(415) 

NSIDES = 2 for l.e.; for the l.e., L must 
be even; For the choice of NE, L, M, 
LLL see suggested values 

2 

LABEL 

(5A10) 

Any alphameric identification 

3 

KEY 10 

(5X,5X, 15) 

Use tolerances preset by code. 

4 

JGEO, JMAT 

(15, 5x, 15) 


5 

TO, TAW, TM, EPS1, TPS2, 
TBG1 , TBG2 

(7F10.5) 

a) If the adiabatic temperature or the 
paint melt temperature vary on the 
surface(s) where h is to be computed, 
leave TAW or TM blank 

b) For l.e.s, EPS1 = EPS2, TBGl - TBG2 

6, 6a, 6b etc. 

TAW(I), . . TAW(K) 

(8F10.5) 

Needed only if the adiabatic wall temp- 
erature vary on the surface where the 
hs are to be computed (more than 1 
card if K > 8). and also if JGEO is not 
2. Omit cards otherwise. 

7,7a,7b,etc. 

TM(1) TM(K) 

(8F10.5) 

Needed only if paint melt temperature 
varies on surface where h is to be com- 
puted. Omit cards otherwise. 

8 

NMAT, RHO, CONDAV, CPAV 

(110, 3F10.5) 

a) Needed only if arbitrary properties are 
inputted (JMAT = -1 ). Omit card 
otherwise. 

b) If material properties vary with temp- 
erature, CONDAV and CPAV are left 
blank. 

c) If material properties are constant with 
Temperature, NMAT is left blank. 

9, 9a, 9b, etc. 

TMATO) TMAT(NMAT) 

(8F10.5) 


10, 10a, 10b, etc. 

CONDT(I), . . . CONDT(NMAT) 

(8F10.5) 

k Cards needed only if properties are temper- 

| ature dependent; omit otherwise. 

1 1,1 1a,1 1b,etc. 

CPT(1), . . . CPT(NMAT) 

(8F10.5) 


12,12a,12b,etc. 

1 3, 1 3a, 1 3b, etc; 

DELX(1), . . . DELX{M) 
DELY(1),. .. (DELY(L) 

(8F10.5) ) 

(8F10.5) ) 

Cards needed only for slab geometry 
(JGEO=0); omit otherwise. 

14 

MCAP, THETA, ALPHA 

( 1 1 0,2 F 10.5) | 

1 

15, 15a. 15b, etc. 

DELX{1), . . . DELX(NWEDGE) 

(8F10.5) 

1 Cards needed only for leading edge geo- 

j metry (JGEO=1 and 2); omit otherwise. 

16, 16a, 16b, etc. 

DELRO ), . . . DELR(L/2) 

(8F10.5) J 


1 7,1 7a, 1 7b, etc. 

X(1), Y( 1 ), . . , X(NN), Y(NN) 

(8F10.5) 

Needed only for arbitrary 2D geometries 
(JGEO= -1 ) ; omit cards otherwise 

18 

JSIDES, KPTS 

(215) 

KPTS is always even 

19, 19a, 19b, etc. 

X(1) X(KPTS) 

(8F10.5) 


20, 20a, 20b, etc. 

T(1 ), . . . T(KPTS) 

(8F10.5) 


21 

RR, QQ 

(2F10.5 

If smoothing not needed use RR = 0.0, 
QQ = 0.001 

22 

EMINF, GAM, PINF, TTOT, 
RGAS.PR 

(6F10.5) 

Card needed only with special option for 
l.e. (JGEO=2); omit otherwise 


TABLE A-2. NOMENCLATURE 


Symbol 

ALPHA 

CONDAV 

CONDT(I) 

CPAV 

CPT{I) 

DELR(I) 

DELX(f) 

DELY(I) 

EMINF 

EPS1 

EPS2 

GAM 

JGEO 

JMAT 

JSIDES 

K 

KEY10 

KPTS 

L 

LLL 

M 

MCAP 

NE 


Effective angle of attack of l.e., positive as in Figure A-2 ('effective" means in a 
cut normal to the model l.e.) 

Average thermal conductivity of model in the temperature range of interest 

Thermal conductivity of model (function of temp.) at the temperature 
TMAT(I), I = 1, NMAT 

Average specific heat of model in the temperature range of interest 

Specific heat of model at temperature TMAT(I) (function of temperature) 

Spatial increments in radial direction for l.e. as indicated in Figure A-2; 

I = 1, L/2 

Spatial increments in x direction of slab geometry as indicated in Figure A-t; 
also for the l.e. as in Figure A-2; I = 1, M 

Spatial increments in y direction of slab geometry as indicated in Figure A-1; 

I = 1, L 


Effective free stream Mach number 

Emissivity of "upper" surface of slab-like geometry or l.e. surface 
Emissivity of "lower" surface of slab-like geometry or of the l.e. surface 
Constant ratio of specific heats of tunnel gas 

Geometry index; selects geometry and, for the l.e.. the type of reduction pro- 
blem, as follows: = -1, arbitrary four-sided geometry; = 0, slab; = 1, l.e. with 
all melt times given; = 2, l.e. with special option with minimum melt time given 

Model material index; selects model material as follows: = -1 arbitrary; = 0, 
Stycast with constant properties, automatically averaged over appropriate 
temperature range; = 2, properties used in NASA submitted problem 

Same as NSIDES 

Number of surface elements on which the hs are to be calculated; = NSIDES 
x M for four sided geometries; 2 x M for l.e. 

Index that selects tolerances: KEY 10= 0, tolerances preset in code 

Number of positions on the model surface(s) at which the times of melt are 
given {< K). Same position number on upper and lower surfaces 

Number of elements through the material (see figures A-1 , A-2, A-3) 

Number of time steps for the time integration of the temperature correction due 
to temperature dependent properties or radiation; if the material properties are 
constant and radiation is neglected, use LLL = 5 

(for slab-like geometry) number of elements along each surface on which the 
hs must be calculated (if 2 surfaces, must be same for both); (for l.e.) number of 
elements on half l.e. surface — see Figure A-2 

Number of elements into which nose of l.e. is subdivided (see Figure A-2); 
must be even 

Number of dominant eigenvalues (substancially less than N) 


Units 

degrees 

BTU/sec ft°F 
BTU/sec ft°F 

BTU/lb °F 
BTU/lb °F 
feet 

feet 

feet 
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TABLE A-2. NOMENCLATURE (CONTINUED) 


Symbol 


Units 

NMAT 

Number of points in thermal properties tables 


NN 

Number of points to specify corners of elements (arbitrary geometry - see 
Figure A-3) 


NSIDES 

Number of sides (1 or 2) in the slab that are heated (a side not heated is 
taken as adiabatic) 


NWEDGE 

Number of elements on the wedge portion of l.e. (see Figure A-2) 


PINF 

Free stream pressure of tunnel gas 

psf 

PR 

(Constant) Prandtl number of tunnel gas 


QQ 

Spline fitting parameter 


RGAS 

Gas constant of tunnel gas 

ft-lb/slug °R 

RHO 

Density of model material (constant with temperature) 

lb/ft 3 

RR 

Spline fitting parameter (= 0 no smoothing; = 1 straight line) 


T,T(1) 

Time(s) of melt; 1 = 1, KPTS 

seconds 

TAW 

Adiabatic wall temperature, constant on the surface(s) where hs are to be 
calculated 

°R 

TAW(I) 

Adiabatic wall temperature, variable or just different constants on upper and 
lower surface; 1 = 1, K 

°R 

TBG1 

Background temperature for radiation from 'upper' surface of slab-like 
geometry or from the l.e. surface 

°R 

TBG2 

Background temperature of radiation from 'lower' surface of slab-like 
geometry or from the l.e. surface 

°R 

THETA 

Wedge half angle of l.e. geometry as indicated in Figure A-2 


TM 

Melting temperature of the paint, if constant on the surface(s) 

°R 

TM(I) 

Melting temperature of the paint, if variable on the surface(s) or even constant 
and different between two surfaces; 1 = 1, K 

°R 

TMAT(I) 

Temperatures in thermal properties table, i.e. temperatures at which a value of 
CONDT(I) and CPT(I) is given; 1 - 1, NMAT 

°R 

TTOT 

Total temperature of tunnel gas 

°R 

TO 

Initial temperature of model 

°R 

X (in cards 19) 

arc lengths at which the times of melt are given (see Figures A-1, A-2, A-3) 

feet 

X (in cards 17) 

X coordinate of corners of elements (arbitrary geometry, see Figure A-3) 

feet 

Y 

Y coordinate of corners of elements (arbitrary geometry, see Figure A-3) 

feet 
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TADIABWALL) npnrprn /?\ /f>> 

MELT TEMPER ) 0RDERED W (!i) 

TIMES OF MELT I ORDERED (a) (?) 

THEIR LOCATIONS x A x B ) KPTS = 


FIGURE A-1 . INDEXES AND INPUTS FOR SLABS 
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EXAMPLE SHOWN: 


M = 5, L = 5, N = 35 

NO. OF CORNERS OF ELEM = (M + 1) (L + 1) = 36 


COORDINATES OF CORNERS OF 



(L+1) (M+1) 


L ELEM * 


FIGURE A-3. INDEXES AND INPUTS FOR ARBITRARY GEOMETRIES 






TABLE A-3. OUTPUT FORMAT OF INPUT DATA AND PRELIMINARY CALCULATIONS 


3 .5*2 THE HEAT TRAnSFfp PRO>M m "*• * OF CT AKGLr#*Fl.E**EMT5 NUMREnEr Cm UMNutSF — 1 0» »l>, <>!» • INV OPT* H*PT B 


12 


TNPuT data 
? «> 
NSIDES L 


l 

NE 


5 

LLL 


N 


?* surface ELEMENTS— 9 BOWS WY 12 COLUMNS 
1 325 ' *13 5*1 6*9 6?* 

\299R 1357* i36«2 1373? 137M 1376^ 

1 *0*9 1*73 1*‘97 1*121 1*1*5 1*169 


OTVES 

698 

1-*0O9 

1*277 


*8 FLEMFnTS 
722 7*6 

13B33 13857 

1 *385 1**93 


NUMBER OF 
EIGENVALUES NE 

n nOMNANT MftoES. • .Rt^UT»E5 
8* * 878 9i>? 12566 

1 708 1 13955 1392* 13953 

1 *6 a ) 1 *6t n 1*710 15697 


15680 W0PD5 OF MEMORY 
1267* 12^82 12890 
11977 1**0 l 1*925 


ECOnOMOM I 7E • • .REOUfF 
LABEL INFORMATION 
0 2 

5*0. 00£o0l*o0* 5gO0P 

toV 12 


CoNiS 1 


DIMENSION OF 5 AND VALUE 
E 1 

n.noroo o. 


M, ' 
i • OOnOO 
IF CONSTANT 
RHO • 14r 


of macros fr™ 
E„ T„ 


?f)000 TOWARDS 1*ft9 5 > 


ftpf'00 


N 

SMI 

CA* 

C*Y 

1 

2-1.3557* 

.00001 

• C 00 77 

2 

1*2. *5072 

000O2 

.0006* 

3 

120.39283 

.00003 

.00055 

* 

120* <19283 

.00003 

.00055 

5 

120.39283 

.00003 

•60055 

6 

120. 39283 

'00003 

.00055 

7 

120.39283 

.00003 

.(006* 

6 

1*2. *5072 

.00002 

.C0077 

9 

2ol.*55T* 

.00001 

o.ooooo 

10 

20 1 .*557* 

.00001 

• 1; 0 0 7 7 

11 

1*2. *5072 

.00002 

.0006* 

12 

120.39283 

.00003 

.00055 

13 

120.39283 

.00003 

.00055 

1* 

120.39283 

,00003 

.00055 

15 

120.39283 

.00003 

.00055 

16 

120.39283 

,00003 

♦ 1006* 

17 

1*2. *5072 

,00002 

.00077 

18 

2*1. *557* 

.00001 

O.CQQOO 

19 

2 .->1.9557* 

.00001 

.10077 

20 

1*2. *5072 

.00002 

.0 1)064 

21 

120.39283 

.0000* 

.C0055 

22 

120.39283 

.0000* 

.00055 

23 

120.39283 

.0000* 

.00055 

2* 

120.39283 

.00004 

.00055 

25 

120.39283 

.0000* 

.0006* 

26 

1*2. *5072 

.00002 

.00077 

27 

21. *557* 

,00001 

o.OOOOO 

28 

284.9nl** 

.00002 

.00038 

29 

2*1. *557* 

.0000* 

.00032 

-><» 

170.26118 

.00005 

*00027 


170.26118 

.00005 

.00027 


' 7o. 26118 

,00005 

,00027 


’Mlfl 

.00005 

,00027 



AAAft* 



. *0 

C.OOOOO 

,00019 


2 ‘BGl, 

*000 5*e.poo’o 5**.oonno 

CP * .22000. 


QUANTITY L : 

-VALUE ACTUALLY NEEDED 
- VALUE USED 


102 

103 

10 * 

105 

106 

107 

108 


28*. 901** 
2*0.78567 
2*0.78567 
2*0.78567 
2*0.78567 
2*0.78567 
20* ,901** 
*02.911*8 


0.00000 
0.00000 
<>.00000 
*.00000 
1.00000 
U.OCOOO 
0.00000 
(J. ooooo 


.00016 
. 0001 * 
. 0001 * 
. 0001 * 
. 0001 * 
. C 00 1 6 
.00019 
C • 0 0000 


.SMOOTHING FACTOR 
FOR UPPER SURFACE 


. 001600 
,00*800 
POSITION * 008O00 


3.6 


6 

(COCO 


AT 
WHICH 

times 

OF 

MELT 

ARE 

GIVEN 


. 010*00 

• 012 r 00 

.013600 

.015200 
.016*00 
.017200 
.01 8nO(J 
.010800 
.019600 


GIVEN 
TIMES 
OF 


3.600000 
•*.600030 

3.600000 

3. 600 000 

MELT 

3.6oO(jJO 
3.609000 

3.600000 
3.60COOO 
3.606000 
3.600030 


1 3.600000 
3.60C 0 00 
TIMES ■*. 600000 
OF 3.60C000 
3.60C0Q0 
3. #<00000 
3.*OCOOO 
3.# 00000 
3.600000 
3.600000 
3.600000 
3.60C 000 


MELT 

USED 


_ SMOOTHING FACTOR 
FOR LOWER SURFACE 


CONSTANT OR AVERAGE 
MATERIAL PROPERTIES: 
v CONDUCTIVITY BTU/SEC-FT-°R 
DENSITY LB/FT3 
SPFrTfTi- heAT BTU/LB-°R 


FOR EACH ELEMENT 
OF GEOMETRY, NUMBERED 
AS IN FIG A-l, A-2, A-3- 


CAX(I) = 


CAT0) = 




V 


0 


N 


tv 


CAX=0 last column 
CAY=0 last row 


SMOOTHING OF EXPERIMENTAL TIMES-OF-MELT 


1st, 2nd, 3rd DERIVATIVES 
{ OF SMOOTHED TIMES OF MELT } 

FP 

O.OoCOQO 
O.CrOOO* 

0 *000000 

0.00000* 

0 .OOOOOn 
C. 0*000* 

0.00009* 

0.0O0000 
0.0*00u* 
o.oroot* 

0.0*002' 

0*000000 


7 

6 

F 

FP 

.001600 

3.6C«0';0 

3.600000 

O.O'COjO 

.00*800 

7.600000 

7.60*000 

0.0*0000 

.008000 

3.600000 

3. *00000 

O.OnOOO* 

.010*00 

3.600000 

3.60C000 

0*0*000* 

♦ 012000 

3.6ii<J0l.O 

3.600000 

0 • 00 0 Ou " 

.013600 

3.6C00L0 

3.690000 

0*00000* 

.015200 

3* 40*000 

7.600000 

0.0*0000 

.016*00 

3.60O000 

3,600000 

O.OoOQu* 

.017200 

3.b0Ot);)C 

3.600000 

o.oroooo 

.01 8')0y 

3.60 OO^C 

3.600000 

0.0*000* 

.0l88 0( j 

3,60*000 

3.60CC00 

0*0*000* 

.019400 

3.600000 

3.600000 

0.0*0000 


FPP 

FPPP 

c , 000000 

O. r p{*nf)A 

o.nooooo 

O.*00n0n 

r . 000000 

n.oOOonn 

r.OOOOQO 

n ."OOOOn 

o.nooooo 

o .foOoo* 

(>.000000 

O .nOpnon 

f. oooooo 

O.onOOO* 

0,000000 

n.noonno 

c.ooocoo 

O.npOOOn 

o.nooooo 

O.rOOOO* 

' .onoooo 

O.nOOPOO 

n.OOCOOO 

O .AftOOOft 

FPP 

FPPP 

c.ooocoo 

o.oooopn 

C.nooooo 

n.rOOooo 

f. 000000 

n.oooonn 

r .000000 

O.nOOOOO 

f .000000 

o . rOOonn 

O. 000000 

O . nooonn 

r. 900000 

rt.nflOonS 

n.nOOCOO 

« .OOOOOn 

( .000000 

fl.nnooo* 

!■ .000009 

o .rooon* 

f-.OOOOoO 

o . n noon* 

f .onoooo 

o>onon* 


UPPER SURFACE 


LOWER SURFACE 


ALL SYMBOIiS AS 
IN FIG. A-l, A-2, A-3 
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TABLE A -4. FORMAT OF OUTPUT OF FINAL RESULTS (SHEET 1 OF 2) 
a) FOR A SLAB 


INTERPRETATION OF PHASt-L H AN6E P*lNT OATA 

GEOMETRY! SLAB HEATEO ON TWO SIDES 

LENGTH « .020 FEET THINNESS * >0r5 FEET 

INITIAL TEMPERATURE s 5*0.00 OEG R 

EMISSIVITY(X) « 0.000 -< y, UPPER SURFACE 

EMlSSIVITY (2) « 0.000 ^ 

BACKGROUND TEMPERATURE (1) ■ 540.00 DEG R * 

BACKGROUND TEMPERATURE ( 2 ) « 5*6.00 OEG R 





upper 

SURFACE 


RATIO -r 

HEAT TRANSFER 







h SIS 

COEFFICIENT 

ELEm 

X 

TIME (PC) 

TEMP (PC) 

TEMP (aw) 

H(SIS) 

F (CORR) 

H 


FEET 

SEC 

UEG R 

DEG K 

B/S*F2*DF 


0/S*F2*DF 

1 

.0016 

3.6000 

719.5000 

1400.000b 

.6671258 

.69871 

.00*9789 

2 

, 00*8 

3.6000 

729.0000 

1400, 000 b 

, 00^5861 

.65056 

.00*9352 

3 

.0080 

3.6000 

767,0000 

1400.000b 

.6k95437 

.59919 

.0057185 

4 

.0104 

3.6000 

840.000c 

1400.000b 

.0139010 

.60592 

.0084230 

5 

. 0 l 2 o 

3.6000 

91*. 0000 

1 * 00 ,0600 

,0194120 

.75985 

.01*7503 

6 

.1136 

3.6000 

992.0000 

1400.000b 

.027020* 

.86085 

.0232606 

7 

.1162 

3.6000 

I 088 .S 0 O 0 

1400.000b 

.0358432 

.89689 

.0321476 

8 

,0164 

3.6000 

1097.5000 

1400.0000 

.0*26141 

.92960 

,03961*0 

9 

.0172 

3.6000 

1117.5000 

1400.0000 

.6*67515 

.92788 

.0433799 

10 

.0180 

3.6000 

1132,5000 

1400 .0000 

.0502244 

.93225 

.0468218 

11 

.0188 

3.6000 

1142,0000 

1400. 000b 

.0526212 

.92816 

.0488411 

12 

.0196 

3.6000 

1 146.5000 

1400.000b 

.0538116 

.92352 

,0496962 




LOWER 

SURFACE 




1 

.0016 

3.6u00 

692.0000 

1400.000b 

.6658467 

.51267 

.0029974 

2 

.0048 

3,6000 

700 • ouOo 

1400.000b 

.0062084 

.47743 

.0029641 

3 

.0080 

3.6000 

726.0000 

1400. 000b 

.0674396 

.42179 

.0031380 

4 

.0104 

3.6000 

762.0000 

1400. 0000 

.6692757 

.31*07 

.0029132 

5 

.0120 

3.6000 

792.0000 

1400.0000 

.0109366 

.26358 

.0028827 

6 

.0136 

3.6000 

825.0000 

1400. 0000 

,0l293?6 

.23035 

,0029790 

7 

.0152 

3.6000 

856,0000 

1400.000b 

.6)49884 

.21489 

.0030710 

8 

.0164 

3.6000 

874.5000 

1400.000b 

.0163084 

.19193 

.0031300 

9 

.0172 

3.6000 

884,5000 

1460 . 0600 

. 01 70 S 44 

.18109 

.0030884 

10 

.OlBu 

3.6000 

892.0000 

1400.000b 

.61 76356 

.17011 

.0030001 

11 

,0l«8 

3.6000 

897.0000 

1*00. OOOU 

.0180279 

.15862 

.0028596 

) 2 

.0196 

3.6000 

900.0000 

1400.000b 

.0182664 

.16677 

.0030*62 


NOTE: X COORDINATE, 

ORDERING OF ELEMENTS, 
SYMBOLS AS IN FIG. A-l. 


94 


TABLE A-4. FORMAT OF OUTPUT OF FINAL RESULTS (SHEET 2 OF 2) 
b) FOR A l.e. SPECIAL PROBLEM 

interpretation of pha st -change paint data 

GEOMETRY i CYLINDRICAL LEADING tUGt FOLLOWED BY wEOGE 

radius r s .00* feet s * .00* feet afterbody length 

theta = lb. 000 degrees alpha * 15.000 

LEES DISTRIBUTION UStD FOR HEaT TRANSFtR COEFFICIENT I THESE LINES APPEAR ONLY 
MACH NUMBER * 8.000 GAMMA » 1,400 | FOR SPECIAL PROBLEM 

INITIAL TEMPERATURE = S40.00 DEG R 

EMlSSlVITV(l) * 0.000 
EMISSIVITY (2) * o.OOO 

BACKGROUND TEMPERATURE(I) = 540.00 DEG R 

BACKGROUND TEMPERATURE (2 ) * 546.00 DEG R 





upper 

SURFACE 


RATIO -jr 

h SIS 

HEAT TRANSFER 
COEFFICIENT 

elem 

X 

TIME (PC) 

TEMP (PC) 

TEMP (Aw) 

R(SIS) 

F (CORR) 

H 


FEET 

SEC 

DEG R 

DEG « 

B/S*F2*DF 


B/'S«F2*DF 

1 

.0005 

*.5Soi 

66U.OOO0 

1284.6702 

.0061361 

.66977 

.00*1090 

2 

.0016 

*.7502 

660.0000 

1265.1762 

. 0061918 

.*8776 

.0030201 

3 

.0026 

5.2005 

660.0000 

1232.4878 

•0?62*?7 

.25981 

.0016219 

4 

.0037 

5.8007 

660.0000 

1205.616V 

.0061907 

.12706 

.0007866 

5 

.00*7 

6. *51o 

660.000P 

1180. *75* 

.0*61*25 

•05206 

.0003198 

6 

.0062 

7.4963 

660.0000 

1168. 7o*y 

. 0*58256 

.03309 

.0001928 

7 

.0072 

7.8821 

660 . O00o 

1168.7041/ 

.0056817 

.03359 

,0001908 

8 

.0082 

8.o959 

660.0000 

1168. 7o4w 

.005606* 

.03370 

.0001889 

9 

.0092 

8.2651 

660*0000 

1168.7040 

.0*55*89 

.*3371 

.0001871 




LOWER 

SURFACE 




l 

,0005 

4.5501 

660.0000 

1296.1852 

.0060290 

.8094 

,00*83*9 

2 

.0016 

4.7502 

660*0000 

1300. OOOu 

.0058669 

.82409 

.00*8349 

3 

.0026 

5.2005 

660 * 0 0 0 o 

1246.1852 

.0056*01 

.85724 

.00*8349 

4 

.0037 

5.8007 

660*0000 

1284.6702 

•0nS*3e8 

.75606 

.00*1098 

5 

.00*7 

6.451 0 

660*0000 

1265.1762 

.Of 531e0 

.56822 

.0030201 

6 

.0062 

7.4963 

660.000c 

12*6.8929 

•Of 56779 

.39646 

.0020132 

7 

.0072 

7.8821 

660*0000 

1246. 892? 

,0"*95*5 

.36683 

.001817* 

8 

.0082 

8.0959 

66() • uOOo 

12*6.8929 

.0**8899 

.3*142 

.0016695 

9 

.0092 

8.2651 

660*0000 

12*6.8929 

.00*8*06 

.32077 

.0015527 


NOTE: ARC LENGTH X, 

ORDERING OF ELEMENTS, 
SYMBOLS AS IN FIG. A -2 
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ADD 40 sec 

FOR COMPILATION 



NUMBER OF ELEMENTS 


FIGURE A-5. COMPUTER TIME REQUIRED BY CAPE 


l.e. SPECIAL PROBLEM 
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3) for the number of dominant eigenvalues NE, use a standard NE = 5. Absolute 
safety is afforded by NE = 10, but the CPU time savings between 5 and 10 are 
some 50% and are therefore attractive. 

4) Since the material properties variations with temperature do not have but a 
tiny effect, great sophistication in the table of the properties is not warranted. 
Use, for example, NMAT = 5. 

5) LLL = 5 has been found adequate. 

6) The tolerances set automatically by the code if KEYIO = 0, are typical, per- 
haps on the tight side, namely at all surface points, simultaneously the h’s differ 
by less thanl% between iterations and the temperature errors, T melt - T ccniputed , 
is less than 1 °R. 

Other points worth of attention are the following: 

7) Each element in the geometry must be attributed the ordering indicated in fig- 
ures A-l , A-2 and A-3; it is in this order that such arrays as DELX, DELY, 
DELR, TAW, TM, X, Y are to be structured. For the slab, one first inputs 
the entire set of (say) TAW etc. for the upper surface, followed by the lower 
surface. Forthel.e., similarly first the ‘upper* side of fig. A-2, sequen- 
tially followed by the lower side. Similarly with the X(KPTS) and T (KPTS) 
arrays. However, for DELX and DELR, both in the slab and in the 1. e. , be- 
cause of symmetry, just one set — say for the upper surface — is to be in- 
putted. When inputting the lower surface of a l.e. , start with the point on the 
axis of symmetry (which therefore is considered a double point). 

8) To do a one dimensional finite-slab geometry or a quasi-one dimensional ar- 
bitrary geometry take M>2 respectively in the slab and the arbitrary geom- 
etry; in other words, ‘double* the given input data in the M direction. 

9) To do a semi-infinite slab or a semi-infinite arbitrary four sides geometry, 
take the ‘lower* surface in fig. A-l and A-3 far enough into the material so 
that there the initial temperature is substantially unchanged up to the maximum 
melt time; then prescribe NSIDE = 1 and the code will treat the ‘lower* surface 
as adiabatic. 

10) In 1. e. cases where the input data are symmetric and a = 0, the code works 
out the entire problem disregarding the symmetry; there is no artifice to take 
advantage of the symmetry. On a slab or a four-sided geometry, to take ad- 
vantage of symmetry of inputs, either in the lateral or in the depth direction, 
just input one-half of the problem. 

11) Appendix B gives the input cards and printouts for two check cases. 
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APPENDIX B 


PROGRAMMER-ORIENTED DOCUMENTATION OF THE CODE 

In this appendix, the details of CAPE organization and structure are assembled. 

A simplified overall flow diagram of CAPE is given in fig. B-l. 

CAPE is organized in a main program and 33 subroutines. The list of the sub- 
routines is given in table B-l together with the function of each subroutine and the sub- 
routine from which it is called. 

The overall logic of the program can now be organized in terms of the subroutines 
as shown in fig. B-2. 

The flow charts for each subroutine are given in sheets B1 - B29. (four routines 
ERFCC, COTAN, ARCOS, and OUT are omitted since they consist of obvious few lines 
of coding). The object deck of each subroutine includes comment cards for each of the 
operations described in the flow chart so that correspondence between the instructions 
and the flow chart can be done at once. 

The set of input cards for two check cases are presented in sheet B-30 while the 
respective outputs are given in table A-4. 

The standard output has been described in Appendix A. For detailed output that 
maps the iterations and eigenvalues and also times the main steps of the calculation, 
the flags LTE and MON must be set equal to six in a statement card at the very begin- 
ning in subroutine PCP SIZE. All the detailed output refers to the iterations performed 
in the subroutine DETRAD and the subroutines that are called from it. 

Self-explanatory diagnostic messages are printed out for the two main failure modes 
(that are intrinsic to the method, rather than to errors in input preparation and in- 
tegrity of the code). One mode is a failure in the inversion if the matrix Gjj is singular. 
The other mode is the failure of the Jennings algorithm to converge, within the maximum 
number of iterations, to the requested number of dominant eigenvalues and eigenvectors. 
The latter mode of failure has never been encountered in normal runs. Presumably the 
corrective action would be to increase the allowed number of iterations, fixed by the 
index NIJ that is set in statement at the outset of subroutine SIZE. The former failure 
mode has been encountered once while exploring the extreme values for NE, the number 
of dominant eigenvalues and eigenvectors. It was found that one slab-like case gener- 
ated a singular G tj matrix for the extreme value of NE equal to 2. The corrective action 
is, of course, to use a reasonable number of NE, from 5 to 10. 
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INPUT: JSIDES, L, M 


0, SLAB 

JGEO } 1, L.E. <t m bs. X) 

2, L.E. (LEES) 

-1, ARBITRARY 
O, STYCAST 

JMAT J 1 , STAINLESS ST. 

2, NASA DEMO. 
-1, ARBITRARY 


T|NIT 

ei,e 2 

T B G, t BG, 2 
T M vs. X (TABLE) 
Taw vs. X (TABLE) 
tM vs. X (TABLE) 

R 




FIGURE B-1 . OVERALL FLOW OF LOGIC IN CAPE (SHEET 1 OF 3) 
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FIGURE B-1. OVERALL FLOW OF LOGIC IN CAPE (SHEET 3 OF 3) 
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TABLE B-1. SUBROUTINES USED IN CAPE 


Name 

Called by 

Main Purpose 

SIZE 

MAIN PROG 

Computes storage locations needed. Compares to number requested. 

PCP 

SIZE 

Reads main input, output initial data. Performs preliminary calculations. 
Controls geometry, material properties and other options. 

SLAB 

PCP 

Computes volumes and conduction shape factors for slab geometry. 

LEAD 

PCP 

Computes volumes and conduction shape factors for leading edge 
geometries. 

ARBIT 

PCP 

Orients individual elements of the arbitrary geometry for computations in 
ELEM. 

ELEM 

ARBIT 

Computes volumes and conduction shape factors between two arbitrary 
quadrilateral elements of the arbitrary geometry. 

SMOOTH 

PCP 

Sets up calculations for DATFT and normalizes the range for the smoothing 
fit. 

SMOFIT 

PCP 

Finds value from table by smoothing-spline interpolation. 

DATFT 

SMOOTH 

Determines Taylor series coefficients for smoothing spline fit. 

LINFIT 

PCP, DETRAD 

Finds value from a table by linear interpolation. 

ATTACK 

PCP 

Finds element number closest to stagnation point and renumbers elements 
as required by LEES. 

LEES 

ATTACK 

Computes ratios h/h S p and T/^VV variation around leading edge for the 
special l.e. problem. 

DETRAD 

SIZE 

Main routines for the inverse problem calculation. Calls eigenvalue and 
matrix routines. Perturbs h to generate influence coefficients. Performs 
iteration on h. Returns the h values. 

IJEN 

DETRAD 

Obtains dominant eigenvectors and eigenvalues of a given matrix (using 
Jennings method, i.e. by simultaneous vector iteration). 

EIGVC 

DETRAD 

Prepares approximate guesses for the eigenvectors to start the Jennings 
algorithm iteration for the zeroth h iteration. 

BFACS 

DETRAD, IJEN 

Factorizes a banded positive-definite matrix. 

BSOLS 

DETRAD, IJEN 

Using the factors of a given banded positive-definite matrix A as generated 
by BFACS solves for X the system AX = Y. 

ORNML 

IJEN 

Carries out the standard Gram-Schmidt orthonormalization of a group of 
vectors. 

HETRAC 

DETRAD 

Sets up coefficient matrix (of conductances) in compact form. 

RVORDR 

IJEN 

Re-orders estimated eigenvalues according to magnitude. 

AORDER 

IJEN, RVORDER 

Sets up permutation indices needed for ordering the eigenvalues. 

DISPLA 

Various 

Prints information, mainly debug special output, in array form. 

LUSOL 

RIMEQF 

Given the results of RDET, substitutes the solution for right hand side. 

PART 

Various 

Prints debug output information and intermediate timing of calculation. 
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TABLE B-1. SUBROUTINES USED IN CAPE (CONTINUED) 


Name 

Called by 

Main Purpose 

RDET 

RIMEQF 

Given a matrix, it decomooses it into two triangular arrays, a lower and an 
upper A = 

RIMEQF 

DETRAD 

Solves system of simultaneous linear equations of type AX = Y for X. 

SCAPRO 

Various 

Computes scalar products of two vectors and adds a value to the result. 

SWITCH 

DISPLA 

Converts columns of a matrix to rows or visa versa. 

SCAPR2 

ORNML 

Computes scalar product of two vectors (without adding a given number 
to the result) 

ERFCC 

PCP 

Computes complementary error function of a given argument. 

COTAN 

ELEM 

Computes the cotangent of an angle. 

ARCOS 

ELEM 

Computes the arc cos of an angle. 

OUT 

SIZE 

Prints final output of problem. 
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START 


PCP 



INPUT: JSIDES, L, M 
JO, SLAB 

JGEO J 1, L.E, (tM vs. X) 
2, L.E. (LEES) 

.-1, ARBITRARY 
O RTYPA^T 

JMAT J 1,STAINLESSST. 

2, NASA DEMO. 
-1, ARBITRARY 


TlNIT 

«i.e 2 

T BG, 1 TbG, 2 
Tm vs. X (TABLE) 
TAW vs. X (TABLE) 
tM vs. X (TABLE) 

R 


INPUT • p 
k vs. T (TABLE)* 
Cp vs. T (TABLE) 


SLAB, 



USE INTERNALLY 
TABUL. PROPS.: p 

k = Vi [k(T | |s] u> + k (T|\/|)] 

Cp = Vi [Cp (T|NIT> + Cp (Tm)1 


1,2 (L.E.) 


■n f^ B,T 


INPUT 

AXi , AX2 , . . ., AXm 
AYj.AYj ay l 


L L 

[PCP 

1 


JJ [— 1 _ 


INPUT 

X,.Y lf X,.Y 2l .. 
(COORDINATES 
OF CORNERS 
OF ELEMENTS) 




H INPUT 

M CAP, ©. « 

Ar, , Ar 2 , . . Ar|_/2 

J . AX,,AX 2 axmwed! 

I— \-~ 




'SINGLE ENTRY IN 
TABLE MEANS CONSTANT 
VALUE 


COMPUTE 


COMPUTE 


COMPUTE 

kjj (pCpV) i 


kjj (pCpV) i 


kjj (pCpV) i 



CONT'D 



l 


2 (LEES) 


FIG. B-2 OVERALL FLOW OF LOGIC ON CAPE IDENTIFIED WITH MOST IMPORTANT 
ROUTINES. (SHEET 1 OF 3) 
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Jl I L_ 







t 


1 


FIGURE B-2. OVERALL FLOW OF LOGIC ON CAPE IDENTIFIED WITH MOST IMPORTANT 
ROUTINES. (SHEET 2 OF 3) 


JL 













t ( 7 ) t 



FIGURE'B-2. OVERALL FLOW OF LOGIC ON CAPE IDENTIFIED WITH MOST IMPORTANT 
ROUTINES. (SHEET 3 OF 3) 
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SHEET B-1 SUBROUTINE SIZE FLOW CHART 











SHEET B-2 (SHEET 1 OF 2) SUBROUTINE PCP FLOW CHART 
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SHEET B-2 (SHEET 2 OF 2) PCP FLOW CHART 









SHEET B-3 SUBROUTINE SLAB FLOW CHARTS 
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SHEET B-4 SUBROUTINE LEAD FLOW CHART 


ENTER 
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SHEET B-5 (SHEET 1 OF 3) SUBROUTINE ARBIT FLOW CHART 



113 



























SHEET B-5 (SHEET 2 OF 3) ARBIT FLOW CHART 
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DO D 















SHEET B-5 (SHEET 3 OF 3) ARBIT FLOW CHART 



CALL ELEM 
FOR A/AX (D_»0 


\ 

f 

RELABEL COORD, 
(e.g. IN ELEM © ) 
7-— A, n — — »B 
©— C, 6 -D 

1 

r 


( CALL ELEM 
l FOR A/AX 0 _*© 


c 


1 + 1 k 


j=j+i k 



CALL ELEM 
FOR A/AX 



f 

RELABEL COORD, 
(e.g. IN ELEM © ) 

8 ■ ■' ■ A, v — — B 

M — — C, 7 ► D 

1 

f 


f CALL ELEM 
l FOR A/AX^_^ © 


1 

r 

(A/AX) 0 



1 

_L + _L 

( A\ 

{ A\ 

0 VAX/0-® 


NO 


NO 




/( CALL ELEM 
l FOR A/AX©_^ X 


1 

f 

RELABEL COORD. 1 

(e.g. IN ELEM (8) ) 

t/YY _» 

o 

II 

p> 

© “B 

rj . - r> 


\ 

f 


CALL ELEM 
V FOR A/AX X _^ ( 


1 

f 

(A/AX) 



1 

"_L + 

1 

(A\ 

(A') 


I V AX>V© 




(A/AX) q 


1 


1 + 
(A) 

1 

(A) 

V AX7 ©^ X 

V AX' 


SET A/AX OF LAST COLUMN 
(e.g. (9) (§) @ (g) ) = 0 
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SHEET B-6 SUBROUTINE ELEM FLOW CHART 
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SHEET B-7 SUBROUTINE SMOOTH FLOWCHART 
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SHEET B-8 SUBROUTINE SMOFIT FLOW CHART 
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SHEET B-9 SUBROUTINE DATFT DESCRIPTION (SHEET 1 OF 2) 


Smoothing Spline Routine 
TITLE: DATFT 

AUTHOR: Antony Jameson, Grumman Aerospace Corporation 

DATE: January 1972 

APPLICABLE COMPUTERS: IBM 1 1 30, 360-370 series, CDC 6600 series 

SOURCE LANGUAGE: IBM 1130 FORTRAN 


PURPOSE: To generate the smoothest possible curve that will pass within specified tolerances of a given set of 

data points. By reducing the tolerances to zero, the curve can be made to pass through the data 
points. By increasing the tolerances sufficiently the smoothest curve becomes a straight-line, least- 
squares fit to the data. 


METHOD: To construct a smoothing spline: 

Let (y j, Xj), i = 1, 2, . . . , n, be the coordinates of the data points. Then construct the curve f(x) 
that minimizes 


I = 



(f" (x)) 2 dx+ R 2 


n 

2 

i=1 


Vi ~ f i 
Q 


2 


where 


fi = f (X;) 


( 1 ) 


and 


Qj = tolerance for i t ^ 1 point 


Here the second derivative f" (x) is used as a measure of curvature. 

The minimizing curve is a piecewise cubic curve that can be expressed over the interval from Xj to 
Xj + -j in the form 

(X — Xj ) 2 (x-X ;) 3 

f(x> = fj + (x — Xj ) f; + — f r + — 6~ f r (2) 

where 

f'.j, fj', and fj" are the first, second, and third derivatives at Xj. 


The coefficients f j, fj', and f'j' are calculated by a method similar to that described by Reinsch 

(Num. Math. 10, 1967, 177-183.). This requires the solution of a set of linear equations with a 
sparse matrix containing five diagonals. The number of computer operations is directly proportional 
to the number of n data points. This is an advantage compared with some other regression techniques 
where the number of the operations may vary as n3. 
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SHEET B-9 SUBROUTINE DATFT DESCRIPTION (SHEET 2 OF 2) 


METHOD: In the limit R-, 0 the minimizing curve reduces to a spline passing through the points. In the limit 

(Continued) R 2 0 the minimizing curve reduces to a straight line giving a weighted least squares fit to the data. 
In the routine R-j and R 2 are determined by a 'smoothing parameter' R through the relations 

R 2 = (t - R ) 4 , Rt = 1 - R2 

so that R = 0 yields a pure spline (zero smoothing) and R -j = 1 gives a straight line (full smoothing). 


The tolerances can be varied from point to point by choice of the weighting parameters Qj. Smaller 
values of Q* L should be used in regions where the curve is expected to have high curvature, to prevent 
the smoothing procedure replacing the true curve by one of larger radius passing outside the data 
points. 

LIMITATION: The curve that minimizes I has zero curvature at the end points: 
f" (x,) - f"(x n ) = 0 

If the true curve is known not to have zero curvature at either or possibly both of the end points, this 
will lead to a systematic error near the end point where the violation occurs. The minimization 
problem, Eq. (1), with free end conditions ought then to be replaced by a minimization problem 
with appropriate constraints at the end points. The existing routine has no provision for this. 

USAGE: CALL DATFT (M, N, S, F, Q, R, A, B, C, D) 

The input data consists of the arrays F ( I ) of values of the dependent variable, S(l) of values of the 
independent variable, and Q(l) of the tolerances to be allowed. S(l) must be a monotone increasing 
or decreasing array. It is also necessary to supply values for the indices M and N, and the smoothing 
parameter R. The routine generates a fit over the interval from S(M) to S(N). R should be a number 
between 0.0 and 1.0; R = 0.0 gives zero smoothing (pure spline); R = 1.0 gives full smoothing 
(straight line least squares fit). 

The routine returns the arrays A(l), B(l), C(l), and D ( 1 ) defined from I = M to I = N. A(l) are the 
fitted values at S(l). B(l), C(l) and D(l) are the first, second, and third derivatives to be used in 
evaluating the cubic curve over the interval from S(l) to S(l+1 ) according to the formula, Eq. (2). 

SUBROUTINES REQUIRED: None. 
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SHEET B-10 SUBROUTINE UNFIT FLOW CHART 
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SHEET B-11 SUBROUTINE ATTACK FLOW CHART 
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SHEET B-12 SUBROUTINE LEES FLOW CHART 















SHEET B-13 SUBROUTINE DETRAD FLOW CHART 


ENTER 


SET UP QUANTITIES 
FOR e=£0&TD PROP. 


SET TO ZERO 
FORCING TERMS IN 
TEMPERATURE, DUE TO 
e * 0 & TD PROP. 


SET FLAGS, # OF 
SURFACE ELEM. INDEX, 
INDECES FOR SPECIAL L.E. 


PRINT TITLE 


CALCULATE INITIAL 
GUESSES hi IF NOT 
SUPPLIED 


CALL HETRAC 

SET UP MATRIX A 
IN COMPACT FORM 


PRINT 


(FOR hj) 

MAJOR ITERATION 


MINOR 

ITERATION 


(FOR Gjj) 


FIRST hj 
ITERATION 



CALL EIGVC 


INITIAL GUESSES 
FOR E & E'S 


CALL RVORDER 
ORDER 

GUESSED E & E'S 


/ WHICH \ 
ITERATION? 


AFTER 
1ST hj 

ITERATION 

CALL IJEN 

OBTAIN 

DOMINANT E & E'S 


L24 











FOR EACH j # 
CALL BOX © 


RETURN 


/MlNOR\ k 
S LOOP 
^^INISHED^X 

jTyes 

CALL RIMEQF ** 
INVERT Gjj MATRIX, 
OBTAIN Ah; 


COMPUTE 


NEW hj 



PICKS MAX 
MELT TIME 


81-84 

CALCULATE 
STEADY STATE 
TEMP T oo 

(NO CALCULATION 
IF TAW UNIFORM) 


CALCULATE 
VECTOR B = 
V T M 1 / 2 (T jnit -T 00 ) 


91-100 

CALCULATE 
SUM =(Tj - T oo) V*M*j 
FOR ELEM. i 


CALCULATE 
FINAL TEMP. 
ERROR FOR 
kth ELEM. 


le SPECIAL PROBLEM * ONLY FOR STAGNATION POINT 
e ' ^ ™UBLtivi ** 1 x 1 |\/) ATR |x DONE EXPLICITLY WITHOUT RIMEQF 
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SHEET B-14 SUBROUTINE IJEN FLOW CHART 
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SHEET B-15 SUBROUTINE EIGVC DESCRIPTION 


EIGVC computes guesses of the eigenvalues, eigenvectors and associated permutation index that are necessary to 
start the iteration in the Jennings method to calculate eigenvalues and eigenvectors. The formulae used for these 
guesses are: 


ith eigenvalue 
ith eigenvector 


R — -i 

A- e +y/(h/2) * in { f <i - i) } 


V T 

vl 


FOR TOP 
HALF 


FOR BOTTOM 
HALF 



h 
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SHEET B-16 SUBROUTINE BFACS DESCRIPTION (SHEET 1 OF 2) 

Given a special, L-banded, positive or negative definite symmetric N-th order matrix, S, decompose it into the 
product: 

S = U T D” 1 U 


where U is an L-banded nonsingular upper triangular matrix with unit diagonal elements, and D is a nonsingular 
diagonal matrix. 

The S matrix is inputted as elements of three one-dimensional arrays. A, B, E. The N elements of A are the main 
diagonal elements of S, the leading (N-L) elements of E are the L-th super- (and sub-) diagonal elements of S. The 
N-1 elements of B are super- (and sub-) diagonal elements of S. *The trailing L/2 elements of E (optionally) define 
main cross diagonal elements of the upper L-th order submatrix of S. In general, later definitions override earlier one, 
e.g. if L = 1 , B (not E) defines the super diagonal elements of S. In the case for which BFACS is intended, most elements 
inside the band are zero. The cross diagonal is installed only if the argument, a, is not zero. The S matrix is topolo- 
gically equivalent to conduction paths in a slab (leading edge if a ^ 0); consider the N = 12, L = 4 example: 


NOSE i ( 

\ ' 


1 ■ 

~ 

■ 9 

2 - 

■ 6 

■10 

3 • 

7 ■ 

•11 

4 ■ 

• 8 ■ 

•12 


Note: 64 = bg = 0 
for the conduction 
problem generally 
but must be explicitly 
made 0 for BFACS 

Note: If nose paths 
are included 1 / e -j 2 
are used (a =£ 0). 



1 

2 

3 

4 

5 

6 

7 

8 9 

10 

11 

12 


— 



, 

1 






— 

1 

a 1 

b 1 

0 

1 e 12 

1 e 1 

0 

0 

0 0 

0 

0 

0 



1“ — 

\ 








2 

b 1 

a 2 1 

e 11 

1 0 

1 0 

e 2 

0 

0 0 

0 

0 

0 



1 — 1 

1 

1 







3 

0 

1 e 11 1 

a 3 

b 3 

,° 

0 

e 3 

0 0 

0 

0 

0 




V 

J_ _ 
0 










4 

e 1 2 

b 3 

a 4 

'© 

0 

0 

e 4 0 

0 

0 

0 

5 

e 1 

0 

® | 3 5 

b 5 

0 

0_ | e 5 

0 

0 

0 

6 

0 

e 2 

0 

0 

' b 5 

a 6 

b 6 

0 1 0 
1 

e 6 

0 

0 

7 

0 

0 

e 3 

0 

1 

1° 

b 6 

a 7 

CT 

O 

0 

e 7 

0 

8 

0 

0 

0 

e 4 

1 0 

0 

b 7 

a 8 ' ® 

0 

0 

e 8 






1 





— -4 — 



_ — 

_ 

9 

0 

0 

0 

0 

® 

0 

0 

(§)| a 9 

b 9 

0 

0 

10 

0 

0 

0 

0 

0 


0 

0 1 b 9 

1 

a 10 

b 10 

0 

11 

0 

0 

0 

0 

0 

0 

e 7 

0 0 
1 

b 10 

a 11 

b 11 

12 

0 

0 

0 

0 

0 

0 

0 

e 8 1 0 

0 

b 11 

a 12 


. 







1 



- 
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SHEET B-16 SUBROUTINE BFACS DESCRIPTION (SHEET 2 OF 2) 

To take advantage of both symmetry and the band form U is stored in a rectangular array of size MID by L, 
where MID > N; the bottom row is used as scratch storage, r^ t . . . , r^, and the unused bottom triangle is 
zeroed out (for convenience in printing only). U appears as: 


STORED ARRAY NON ZERO ELEMENTS OF U MATRIX 



i 

2 

3 

4 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 


— 



— 1 


H 











— 1 

1 

u n 

u 12 

u 13 

u 14 

1 

1 

U 11 

u 12 

u 13 

u 14 








2 

U 21 

u 22 

u 23 

u 24 

2 


1 

U 21 

u 22 

u 23 

u 24 







3 

U 31 

u 32 

u 33 

u 34 

3 



1 

U 31 

u 32 

u 33 

u 34 






4 

U 41 

u 42 

u 43 

u 44 

4 




1 

U 41 

u 42 

u 43 

u 44 





5 

U 51 

u 52 

u 53 

u 54 

5 





1 

U 51 

u 52 

u 53 

u 54 




6 

U 61 

u 62 

u 63 

u 64 

6 






1 

U 61 

u 62 

u 63 

u 64 



7 

U 71 

u 72 

u 73 

u 74 

7 







1 

U 71 

u 72 

u 73 

u 74 


8 

U 81 

u 82 

u 83 

u 84 ' 

8 ' 








1 

U 81 

u 82 

u 83 

u 84 

9 

U 91 

u 92 

u 93 

0 1 

9 ! 









1 

U 91 

u 92 

u 93 

10 

u 10,1 

u 10,2 

0 

0 

10 | 










1 

u 10,1 

u 10,2 

11 

u 11,1 

0 

0 

0 

11 











1 

u 11,1 

12 

r 1 

r 2 

r 3 

r 4 

12 












1 


The N elements of D“^ are stored in an N-array. For the usual case of S being either positive-or negative-definite, 
these elements are all positive or all negative, respectively. However the routine will "work” provided only that the 
leading N principle minors are non-zero. For details see the following article which guarantees high accuracy only for 
the definite cases of usual interest: "Symmetric Decomposition of Positive Definite Band Matrices", R.S. Martin, 

J.H. Wilkinson, C. 1/4, LINEAR ALGEBRA - HANDBOOK FOR AUTOMATIC COMPUTATION, VOLUME II, 
Springer-Verlag, 1971. 
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SHEET B-17 SUBROUTINE BSOLS DESCRIPTION 


Given the product form decomposition of an L-banded symmetric matrix, S = uTq' 1 U, as calculated by the 
BFACS subroutine, BSOLS solves a system of N linear equations with M right hand sides: 

S {Yj Y 2 - • -YM} = { Y, Y 2 - ■ YM| 

The routine sumply carries out the standard forward substitution phase: 
z = U-T y 

followed by the standard backward substitution phase: 
x = U’ 1 0 z 

The only unusual aspect is the rather unorthodox storage scheme which is described in the documentation for 
subroutine BFACS. This scheme is necessary to exploit the banded symmetric form of S in the most efficient way 
in terms of computer memory. For details see: R.S. Martin, J.H. Wilkinson, "Symmetric Decomposition of Positive 
Definite Bank Matrices", in: Linear Algebra- Hand book for Automatic Computation, Volume 1 1 , C. 1/4, Springer- 
Verlug, 1971 
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SHEET B-18 SUBROUTINE ORNML 


Classical Gram-Schmidt orthonormalization of a set of M linearly independent vectors y! ,"y 2 
in the order Lj , L 2 , L 3 , . . L|yj 











SHEET B-19 SUBROUTINE HETRAC FLOW CHART 
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SHEET B-20 SUBROUTINE RVORDER 


QUANTITY 

SYMBOL 

INPUT/OUTPUT 

DIMENSION 

EIGENVALUES 

R 

IN + OUT 

R(MM) 

EIGENVECTORS 

V 

IN + OUT 

V(MID, MM) 

RANK VECTOR 

K 

OUT 

K(MM) 

PERMUTATION VECTOR 

L 

OUT 

L(MM) 

EIGENVECTOR DIMENSION 

N 

IN 

_ 

NUMBER OF EIGENVECTORS 

MM 

IN 


DIMENSION OF ARRAY 

USED TO STORE EIGENVECTORS 

MID 

IN 




I 
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SHEET B-21 SUBROUTINE AORDER DESCRIPTION 

PURPOSE: ORDER A SET OF REAL NUMBERS 
CALLING SEQUENCE: 

CALL AORDER (A, N, IPERM) 



NAME 

DIMENSION 

DESCRIPTION 

INPUT 

A 

A(N) 

ELEMENTS TO BE ORDERED 


N 


NUMBER OF ELEMENTS = INI 
N>0 INCREASING ORDER 
N < 0 DECREASING ORDER 

OUTPUT 

IPERM 

IPERM(N) 

ORDER VECTOR - 


SPECIFIES THE SEQUENCE 
OF ELEMENT INDEX NUMBERS 
WHICH WILL PRESENT 
A AS AN ORDERED SET, 
i.e. 

DO 100 1 = 1, N 
100‘WRITE (6,1) A (IPERM(I)) 

1 FORMAT (F 10.5) 

Wl LL LIST A AS AN ORDERED ARRAY 


AORDER CALLS NO OTHER SUBROUTINES 
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SHEET B-22 SUBROUTINE DISPLA DESCRIPTION 


TITLE: DISPLA — Prints scalars, vectors, rectangular matrices, packed symmetric matrices, and Hessenberg 

matrices. 

AUTHOR: M.J. Rossi 

DATE: September 1973 

APPLICABLE COMPUTERS: IBM 360/370; CDC 6000 SERI ES 

SOURCE LANGUAGE: FORTRAN IV 


PURPOSE: To simplify printing of mathematical types of data structures in an easily read format which allows 

titles and index labels. 


METHOD: FORTRAN looping and write statements which indexes and addresses arrays according to their type. 

USAGE: Call DISPLA (X, NFILE, TITLE, KAR, KIND, NROWS, NCOLS, MID). 


X — Input — Array of one or more values to be printed 
NFI LE — Input — FORTRAN unit for printing. 

TITLE Input ■ Vector of KA® ch 3 rectors used es + '*!s. 

KAR - Input - Number of characters in above string. 

KIND - Input - Type of mathematical data structure: 

= O scalar (or vector printed on one line with no index) 

= 1 vector of |NROWS| elements, indexed 

= 2 Rectangular IN ROWS! by NCOLS matrix — Dimension (MID. *) 
= 3 Packed Symmetric matrix of order | NROWS I 


1 2 4 

2 3 5 
4 5 6 


— lower triangular partial rows if NROWS positive 


1 2 3 

2 4 5 

3 5 6 


— lower triangular partial columns if NROWS negative 
4 — Transposed Hessenberg matrix of order NROWS — Dimension (MID, MID) 


NROWS 

— Input 

NCOLS 

— Input 

MID 

- Input 


— Number of elements if KIND = 0 or 1 

— Number of rows if KIND = 2 

— Matrix order if KIND = 3 or 4 

— Number of columns if KIND = 2 

— Ignored otherwise 

— Matrix Dimension if KIND = 2 or 4 

— Ignored otherwise 


SUBROUTINE REQUIRED: 


SWITCH 



SHEET B-23 SUBROUTINE LUSOL DESCRIPTION 


Given the factorized product form of A = P * L # U, LUSOL solves a linear system A x = y. The solution vector is 
obtained in two steps: (1) z = (L _1 (pT y)) f and then (2) x = U _1 z. A good discussion of the details may be found in 
the following article: “Solution of Real and Complex Systems of Linear Equations/' H. J. Bowdler, R. S. Martin, 

G. Peters, and J. H. Wilkinson, C. 1/7, pp. 93-110, Linear Algebra - Handbook for Automatic Computation, 

Volume II, Springer-Verlag, 1971 . 
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SHEET B-24 SUBROUTINE PART DESCRIPTION 


TITLE: PART — Prints standard 1 20 character labels at the top of the next page. 

AUTHOR: M.J. Rossi 

DATE: September 1973 

APPLICABLE COMPUTERS: IBM 360/370; CDC 6000SERIES 

SOURCE LANGUAGE: FORTRAN IV with 2 Assembler Language Subordinate Subroutines. 

PURPOSE: To make it convenient to produce standard printed labels with "part” numbers, date, and running CPU 

time on the top of the next page. Also, prints short line on next line with just CPU time for intermediate 
timing. 

METHOD: On the first printing entry for a given computer run the Date Subroutine is invoked and an 8*character 

field of an internal word is stored with the date in the form: "KK/LL/MM", where KK is the index 
number for the month, LL is the day of the month, and MM is the last 2 digits of the year, e.g., 

March 15,1973 3/15/73. Also, at this time, the SECOND subroutine is invoked to both 

establish the zero time point and to set the units to hundreds of a second. Then the first printed page 
heading is given with a zero time and a PART number 1 reported. Subsequent printing calls will give 
the time as: NN.I I .JJ where NN is the number of minutes elapsed, 1 1 is the number of seconds, and 
JJ is the number of hundredths of seconds. The PART number is incremented by one for each 
printing call. There are two fields of alphameric information for the full printing mode which are under 
control of the user: (1) The first is a 40 character LABEL field which is set upon calling PART in the 
non printing mode, (2) The second is a 48 character field which is supplied on a full printing call. 

There is also a partial printing mode which simply results in the appearance on the next line of an 8 
character field of user supplied TITLE along with running CPU time. 

USAGE: Call PART ('XX. . .X'. I \ 

'XX. . .X’ — Input = Alphameric string of either 8, 40, or 48 characters depending on the value of L. 

L — Input — FORTRAN unit for printing, if positive 

— If zero, simply sets 40 character LABEL field and returns 

- If negative, prints 8 character TITLE — 'XX. . .X' — and CPU time on next line and 
increments PART number. 

- If positive, prints DATE, TIME, 40 character LABEL, 48 character TITLE, Part Number 
and spacers with standard notation. 

SUBROUTINES REQUIRED: DATE, SECOND 


t 
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SHEET B-25 SUBROUTINE RDET DESCRIPTION 


Given a square matrix. A, stored in a Fortran double array, decompose it into the product: 

A = P * L * U 

where P is a permutation matrix, L is lower triangular and U is upper triangular. The algorithm includes implicit row 
scaling and partial pivoting while providing a test for singularity of A. For details, see reference given for subroutine 
LUSOL. 
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SHEET B-26 SUBROUTINE RIMEQF 


Solves the system of linear equations, A x = y. The first step involves decomposing A as the product: A = P # L # U 
which is done by a call to subroutine RDET. This followed by the standard method of forward substitution and 
then backward substitution performed by a call to subroutine LUSOL. This is - H uivalent to solving a pair of 
triangular systems using L and then U. The total procedure may be viewed as a variant of Gaussian elimination as' 
described in more detail in the reference given for subroutine LUSOL. 
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SHEET B-27 SUBROUTINE SCAPRO 

SCAPRO adds to a quantity the inner product of two vectors X, Y stored as equally spaced words in Fortran arrays. 

L 

SUM = SUM + 2 X (IX . (J-1) + 1) * Y (IY) . (J-1) + 1) 

J=1 
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SHEET B-28 SUBROUTINE SWITCH DESCRIPTION 


PROGRAM TITLE: Utility routine for re-arrangement of certain triangular arrays 


SUBROUTINE NAME: SWITCH INDEX: 

ANALYST: F. Nolan 

PROGRAMMER: F. Nolan DATE: 

DOCUMENTATION AUTHOR: F. Nolan DATE: 

SOURCE LANGUAGE: FORTRAN IV 

APPLICABLE COMPUTERS: IBM Systems 360, 370; CDC 6000 series 

REVISION: DATE: 


12 . 6 . 0.1 

June 15,1967 
June 20, 1973 


PURPOSE: To provide a convenient conversion between two common arrangements for the storage of triangular 

(and symmetric) matrices. 

ANALYTIC DESCRIPTION: The routine makes systematic use of transpositions, i.e., interchanges of two array 

elements. It is a well known result in permutation theory that every permutation can be represented 
as a product of transpositions. 


PROGRAM DESCRIPTION: There is no loss of generality in assuming that the input matrix is of lower triangular 

form. It is natural to store such matrices by row or by column. Both arrangements are illustrated for 
s matrix of order 5. The 'understanding is thet the 3} e!smsr' + ovampio ^eeinnori nncitmn Q 


using row storage, and position 1 1 using column storage. 
Row Storage 
1 

2 3 

4 5 6 

7 8 9 10 

11 12 13 14 15 


Column Storage 
1 

2 6 

3 7 10 

4 8 11 13 

5 9 12 14 15 


Given a lower triangular or symmetric matrix, stored in either fashion, SWITCH can re-arrange it to the 
other form. The re-arrangement is carried out "in place" in the sense that no auxiliary array is required. 
For an input matrix of order m, the transformation is performed in approximately !4m 2 transpositions. 
There are no rounding errors. 


PROGRAM RESTRICTIONS: 

INPUT PARAMETERS: 

FORTRAN Name 

A 

M 


The matrix must be of order at least 3. 

Description 

Singly-dimensioned real array containing the matrix to be re-arranged. 

Order of matrix A. If M is given positive, conversion is from row to column 
storage, if M is given negative, conversion is from column to row storage. 


OUTPUT PARAMETERS: 

A Matrix in re-arranged order. 

CALLING SEQUENCE: 


CALL SWITCH (A, M) 
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SHEET B-29 SUBROUTINE SCAPR2 


SCAPR2 CALCULATES THE INNER 

PRODUCT OF TWO VECTORS STORED 

AS EQUALLY SPACED WORDS IN FORTRAN ARRAYS. 
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SHEET B-30 INPUT DATA LISTING FOR TWO CHECK CASES 


OOOOOOOOOllllllllll 2222222222333333333 34 44 4 444 44455555555556 6666 66666 7 7777777777 
012 34 5678901234 56 78901234 5678 901 23456 78 90 1234 56789012 34 56 78 90 12 34 56 78 90 1234 56 7 89 

2 10 9 10 5 

LABEL INFORMATION LE LEES ALPHA# 1 5 CP 
- 1 0 0 

2 0 1 


540.0 

o 

• 

o 

660.0 

0.0 

o 

• 

o 

54 0.0 

540 • 0 



10 15.0 

15.0 






0.00 1 

0.001 

0.001 

0.001 





0.0004 

0.0008 

0.0008 

0 . 0008 

0.0012 




2 

18 







0 .000523 

0.001570 

0.002617 

0.003664 

0.0047 1 1 

0.005735 

0 .006735 

0.007733 

0.008733 

0 . 000523 

0.001570 

0.002617 

0.003664 

0.0047 1 1 

0.005735 

0 .006735 

0 . 007733 

0.008733 







4.55 

4.75 

5.20 

5.80 

6.45 

7.20 

7.72 

8.0 0 

8.18 

4.55 

4.75 

5.20 

5.80 

6.45 

7.20 

7.72 

^.00 

8.18 







0.0 

0.001 







8.0 

1 • 4 

125.0 

1300.0 

53.53045 

0.72 




0000000001 1 1 1 1 1 111 122222222223333333333444444444455555555556666666666777777777 

0 123456789012345678901 23456789012345678901 23456789 

0 1 234567890 1 234667890 ! 234567 

2 

9 1 2 

1 0 5 






LABEL information 

SLAB TWO 

SIOES CP 





-1 

0 0 







0 

0 2 







54 0.0 

1400.0 

0.0 

O 

. 

o 

o 

» 

o 

540.0 

540.0 


719.5 

729.0 

767. 0 

840.0 

914.0 

992 • 0 

1058.5 

1 097.5 

1117.5 

1132.6 

1 142. 0 

1146.5 

692.0 

700.0 

726.0 

762.0 

792.0 

825.0 

856.0 

874 .5 

884.5 

892 . 0 

897.0 

900.0 

0.0032 

0.0032 

0.0032 

0.0016 

0.0016 

0.0016 

0.0016 

0.0008 

0 • 0008 

0.0008 

0 . 0008 

0.0008 





0 .00025 
0.00025 

0 .0005 

0.0007 

0.0007 

0.0007 

0.0007 

0.0007 

0.0005 

2 

24 







0.0016 

0.0048 

0 • 0080 

0.0104 

0.0120 

0.0136 

0.0152 

0.0164 

0.0172 

0.0180 

0.0188 

0.0196 

0.0016 

0.0048 

0.0080 

0.0104 

0.0120 

0. 01 36 

0.0152 

0.0164 

0.0172 

0.0180 

0.0188 

0.0196 

3.6 

3 • r> 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 

5.6 

3.6 

3.6 

3.6 

3.6 

3.6 

3.6 


0.0 0.001 


U.S. GOVERNMENT PRINTING OFFICE: 1974- 635*042/6 
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